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Abstract 


Panel flutter is a form of dynamic aeroelastic instability resulting from the in- 
teraction between motion of a aircraft structural panel and the aerodynamic loads 
exerted on that panel by air flowing past one of the faces. It differs from lifting- 
surface flutter in the sense that it is not usually catastrophic, the panel’s motion 
being limited by nonlinear membrane stresses produced by the transverse displace- 
ment. Above some critical airflow condition, the linear instability grows to a limit 
cycle. 

The piesent investigation studies panel flutter in an aerodynamic regime known 
as ‘free molecule flow”, wherein intermolecular collisions can be neglected and loads 
are caused by interactions between individual molecules and the bounding surface. 
After collision with the panel, molecules may be reflected specularly or reemitted in 
diffuse fashion. Two parameters characterize this process: the “momentum accom- 
modation coefficient”, which is the fraction of the specularly reflected molecules; and 
the ratio between the panel temperature and that of the free airstream. This model 
is iclevant to the case of hypersonic flight vehicles traveling at very high altitudes 
and especially for panels oriented parallel to the airstream or in the vehicle’s lee. 
Under these conditions the aerodynamic shear stress turns out to be considerably 
larger than the surface pressures, and shear effects must be included in the model. 
This is accomplished by means of distributed longitudinal and bending loads. The 
foimei can cause the panel to buckle. In the example of a simply-supported panel, 
it turns out that the second mode of free vibration tends to dominate the flutter 
solution, which is carried out by a Galerkin analysis. 

Several parametric studies are presented. They include the effects of 1) temper- 
ature ratio; 2) momentum accommodation coefficient; 3) spring parameters, which 
are associated with how the panel is connected to adjacent structures; 4) a. param- 
eter which relates compressive end load to its value which would cause classical 
column buckling; 5) a parameter proportional to the pressure differential between 
the front and back faces; and 6) initial curvature. The research is completed by an 
investigation into the possibility of accounting for molecular collisions, which proves 
to lie infeasible given the speeds of current mainframe supercomputers. 
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Chapter 1 


INTRODUCTION 


The panel flutter phenomenon is a form of dynamic instability resulting from the 
interaction of the motion of a panel with the aerodynamic loads generated by this 
motion. It differs from classical flutter in the sense that it is not usually catastrophic, 
the motion of the panel being limited by a nonlinear membrane stress induced by 
the transverse displacement. In general the flutter motion corresponds to a limit 
cycle, which is not necessarily simple harmonic. This subject was intensely studied 
during the 60’s and 70’s, and there are several works published in the area. The 
definitive reference seems to be the 1975 monograph by Dowell while a more 
recent survey by Reed, Hanson and Alford ® addresses the requirements directly 
related to the National Aerospace Plane. Several other comprehensive historical 
surveys of the liturature are available^ - ^, including Reference [9] which treats the 
applications of the methods and results in design practice. 

From a historical point of view, the first mention of panel flutter appears to lx* 
related with failures that had occurred on the German \ 2 missile* during World 
War II ^ The problem was initially regarded as an interesting but not funda- 
mentally different aspect of the general field of aeroelasticity. However, the use* of 
the standard methods for aeroelastic analysis (mainly linear methods at the time) 
usually resulted in substantial disagreement between theoretical and experimental 
data. This led to many studies, which disclosed that one of the main difficulties was 
associated with the rather inaccurate techniques and uncontrolled test, conditions 
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(by present standards) of early experiments. On the other hand, from the theoret- 
ical point of view, for some time it was not clear if the fluttering motion should be 
modeled by “standing” waves or “traveling” waves. This may all have started from 
the experience that flags, as well as skin panels made of thin paper, flutter at quite 
low speeds and the motion resembles traveling waves flQl . 

These early difficulties have been overcome, and it has been demonstrated that 
the data from theory and experiment are in satisfactory agreement for a wide range 
of parameters. Among the aspects inherent to panel flutter which are not ade- 
quately described by the usual linear, inviscid aeroelastic models are the structural 
nonlinearities. As the panel bends it also stretches thereby inducing a tension in 
the panel. As mentioned before, the flutter motion is a limit cycle, representing a 
balance between the (unstable) linear panel and airloads and that induced tension, 
which increases the effective panel stiffness. Viscous shear flow effects can also be 
of importance, one of the reasons being that the boundary layer region next to the 
panel can significantly modify the effective Mach number and dynamic pressure 
of the flow. Moreover, the shear stress itself may drastically modify the stability 
boundary if it is larger or comparable to the pressure, as in the case of hypersonic 
flow in a rarefied atmosphere. 


1.1 Physical Nature of the Problem 

There are several types of aeroelastic instabilities which may occur, and these may 
be classified as static (“divergence”) or dynamic (“flutter”) instabilities. In the 
case of panel flutter, divergence is generally associated with subsonic flow l 11 !, while 
flutter usually only happens for supersonic flows. However, it is possible for flutter 
to occur in the subsonic regime in the case of a long panel resting on a contin- 
uous elastic foundation^). i n fact, in this case a divergence-like instability sets 
in first, but it is a very mild one and the panel response only becomes significant 
when the airspeed approaches that for a flutter-like instability. The much larger 
number of experimental studies in the supersonic flow regime [1 3)-[19] reflects this 
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fact that divergence is normally of less technological importance than flutter. Fi- 
nally, it is worth pointing out that for flat panels, flutter in the low supersonic flow 
regime often corresponds to a single-degree-of-freedom instability, while for Mach 
numbers approximately larger than y/2 flutter results from coupling between modes 
(“frequency coalescence”). In the latter case it is clear that any parameter that sig- 
nificantly changes the panel modes and/or frequencies may have a significant effect 
on the stability boundaries. 


To better understand the physical nature of the problem, consider a wind-tunnel 
experiment for a flat panel at supersonic speeds, with no applied in-planc loads. 
Once the flow is established at some specified Mach number, the flutter boundary 
is searched for by increasing the fluid dynamic pressure q (the Mach number is 
held at the initial value). For small values of q below the flutter boundary, random 
oscillations may be observed with dominant frequency components near the lower 
panel natural frequencies: the panel is responding to pressure fluctuations in the 
turbulent boundary layer, but it acts as a mechanical filter and responds piimaiily 
to those frequencies in the pressure field near the structural natural frequencies. The 
maximum oscillation amplitude is very small. When the flutter boundary is reached 
at some critical dynamic pressure, q cr , the oscillation becomes essentially peiiodit, 
with the amplitude increasing with any further increase in the dynamic pressure. 
Actually, the experimental value of q CT is more a matter of definition than a result 
obtained from some precise measurement. An arbitrary but reasonable convention 
is to define the stability boundary as a flutter band in dynamic pressure ' '• The 

upper limit can then be taken as the lowest value of q at which a sustained oscillation 
was observed with an amplitude magnification of 3 to 5 times over that observed 
at lower values of dynamic pressure. On the other hand, the lower limit may lie 
assumed to correspond to the highest value of q for which no oscillations with regular 
frequency or with a significant increase in response occurred. The typical etioi in 
q cr is of the order of 10%. 


There are very few works which tried to investigate failure mechanisms, but at 
least two such mechanisms can be identified and have occurred in practice. The fiist 
one has to do with the stress amplitude due to flutter exceeding the yield stress of the 
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panel material over a major portion of the structure. In this case the failure can be 
catastrophic and very rapid. On the other hand, small stress levels can be associated 
with fatigue failure. Actually, if the stress amplitude and frequency of oscillation 
are known, one may use a conventional fatigue curve (stress vs. number of cycles 
to failure) to estimate the fatigue life for the fluttering panel [20,-1]^ Interesting 
enough, in a study by Dowell [20 1 it is shown that whereas a thinner panel will 
flutter at a lower dynamic pressure than a thicker one, it will be able to exceed 
its flutter dynamic pressure by a greater amount (on a percentage basis) than a 
thicker one for a given fatigue life. Another conclusion from the same work is that 
a significant reduction in panel thickness may be possible if the panel is designed 
for finite rather than infinite (no flutter) fatigue life. Results obtained by Xue and 
Mei* 2 ^ agree with this finding. 


1.2 Methods of Analysis 

There are very few exact solutions, as it is always the case in any area, and most 
of them use the method originally employed by Hedgepeth [22,23] This method 
of solution is applicable to isotropic panels, as in the works by Hedgepeth, or to 
sandwich/orthotropic panel S t 6 U 24 H26] Naturally, these analyses are valuable to 
compare approximate solutions against, but they apply only to the linear problem 
and were developed only by assuming high supersonic flows. 

The most usual approximate method of analysis has been to apply Galerkin ’s 
method (271 to the governing partial differential equation(s) of motion. It is worth 
mentioning that the Galerkin method and the Rayleigh-Ritz method are equivalent 
in most cases. The main differences are that the latter method involves the station- 
arity of the total energy of the system, and that the so-called trial functions only 
need to satisfy the forced boundary conditions, i.e., those conditions directly related 
to the edges’ displacements (deflection and slope). Also, some doubt was initially 
cast on the use of the Galerkin method since spurious flutter results were obtained 
for membrane flutter at high supersonic speeds when exact analyses showed no 
flutter. This problem was studied by Ellen t 28 ^, who showed that the instability is 
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always associated with the coalescence of the two highest eigenvalues. Since only the 
lower eigenvalues have been determined with any accuracy, it becomes clear that at 
any stage more terms must be taken to check the instability. Therefore an infinite 
number of modes is in fact needed. This leads to the conclusion that, although a 
finite value is found for q CT for a given analysis, it has no meaning because of the 
manner in which the instability occurs. 

Once Galerkin’s method is applied to the equation of motion, a system of ordi- 
nary differential equations is obtained, which can be solved in several ways: 

Characteristic equation: when the system of equations is linear, and in the ab- 
sence of any structural damping and nonhomogeneous terms as pressurization 
loads, it is possible to write: [29]— f40] 

|M] {$} + [*]{,} + A [4 <,} = <0> , (1) 

where [M], [A], and [A] are the inertia, stiffness and aerodynamic matrices, 
respectively, and A is a nondimensional parameter proportional to the dynamic 
pressure. Note that the aerodynamic matrix [A] may depend on Mach number, 
as well as other parameters, depending on the aerodynamic theory being used. 
The generalized displacements { 5 } are then assumed to be of the form {<f>} e iu,t , 
where the frequency u is in general complex, i.e., u = u; R + iuj j. After this 
substitution, a nontrivial solution of equation ( 1 ) is obtained by setting the 
determinant of the resultant coefficient matrix of {<f>} equal to zero. This 
yields the characteristic equation of the problem, which can be solved for u> 
if the parameter A is specified. For A = 0 the values of to can be shown to 
correspond to the natural frequencies of the panel in a vacuum. Increasing A 
changes the frequencies smoothly until a critical value (A cr ) is reached, when 
two frequencies coalesce. Any further increase in A makes these frequencies 
complex conjugate pairs, indicating that at least one of the modes of oscillation 
is unstable, and thus defining the flutter boundary. Actually, this problem can 
be looked upon as an eigenvalue problem and solved as such. This approach 
seems to be numerically more efficient than to solve the characteristic equation 
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itself if the total number of assumed modes is large. Finally, a traveling wave 
solution can be obtained in a similar way, as in the papers by Dowell [41,42) 

Harmonic balance method: this method is used when the system of equations 
is nonlinear. It is based on the assumption that the flutter motion consists 
of a steady periodic vibration, the fundamental harmonic being predominant, 
thatis,t 43 H46] 

q n (t ) = a n smut + b n cos ut . (2) 

Substituting this into the system of equations and balancing terms of the 
first harmonic (terms multiplying sinwf and cos ut), one obtains 2 N algebraic 
equations for the 2 N + 2 unknowns {X } = {dx , iq , a 2 . • • , b^} T plus A 

and u. It should be recalled that neither A cr nor the relation between A and 
u are known beforehand. One may replace the first two coefficients ai and 
bi in {X} by A and u. Furthermore, since the phase angle is not important 
in a steady-state solution, one can set ax = 0 and bi equal to a certain small 
number, which will set the general level of the vibration amplitude corre- 
sponding to given values of A and u. This reduces the system the equations 
to a determinate form. It should be said that in the recent work by Yuen 
and Lau^ the so-called incremental harmonic balance method is used, with 
higher harmonics being used to carry out the analyses. 

The stability of the solution obtained may be studied by giving a small per- 
turbation to the limit-cycle solution as follows: 

q n {t) = [a n + £„(<)] sin ut + [b n + cos ut . (3) 

Substituting this into the system of equations, and then studying the behavior 
of f n (f) and T/ n (f), one can find out the stability of the limit-cycle solution. 
Note that both variables are small perturbed functions, so that one can neglect 
higher order terms and thus deal only with linear differential equations in £ n (f ) 
and rj n (t). 
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Perturbation method: this method obviously applies to nonlinear cases, the so- 
lution being studied for values of A in the vicinity of the linear value A cr . A 
very clear application of this method can be found in the work by Morino ^8], 

Time integration: in this case the generalized coordinates are obtained by direct 
numerical integration [45] t [49]-[61]^ a p ow i n g the study of nonsimple harmonic 
limit cycle oscillations or even chaotic motion ^1. p or a given initial distur- 
bance, the panel motion may decrease with time (if the motion is stable in a 
linear sense) or may increase with time until a nonlinear but “stable” limit 
cycle is reached. One drawback of this approach is that if all damping is 
eliminated from the system the limit cycle cannot be reached. This is a result 
of treating the problem in the initial value formulation rather than seeking, 
a priori , sustained oscillatory solutions. With regard to convergence, pres- 
ence experience suggests that at least six modes must be used. However, if 
in-plane or static pressure loadings induce a significantly large tension in the 
panel, even more modes may be required. The number of modes required for a 
given desired accuracy also increases if the panel aspect ratio (length/ width) 
increases. 

Another method of analysis which became very common in the last two decades 
corresponds to the finite element method ^1, recognition of the first application 
usually being given to Olson The basis of this method is to divide the panel in 
several discrete elements, the displacement(s) at the nodal points (i.e., points where 
different elements touch each other) being the unknown variables. One of the biggest 
advantages of this method is that it can be applied to problems with practically any 
geometrical boundary conditions on any or all the sides of the panel once the ele- 
ments are defined. In the case of Galerkin’s method a different set of functions must 
be chosen to represent the deformations which correspond to each set of boundary 
conditions. One issue regarding convergence in the case of plate elements has to do 
with the quantities that are taken to be continuous across the interfaces between 
elements. If the interface compatibility conditions are exactly satisfied the element 
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is said to be “conforming”, this being a desired property. A “non-conforming” el- 
ement usually implies that normal slopes are discontinuous across the interfaces, 
while “hyper-conforming” elements use components of the curvature tensor as gen- 
eralized coordinates. In the latter case a degree of continuity that is not required by 
compatibility is introduced, slowing down the convergence properties if compared 
to a solution using the same mesh size and the same degree of polynomial approxi- 
mation but strictly conforming 169 '. Additionally, hyper-conforming elements raise 
difficulties when elements of different thicknesses or Young moduli meet at a ver- 
tex. Two other points that concern only the nonlinear formulation are discussed 
by Sarma and Varadan' 64 '. The first one is related to the inclusion of the in-plane 
displacement term in the strain-displacement relation, since the order of this term 
is the same as the (w, x ) 2 term. Here w, x represents the derivative of the transverse 
displacement with respect to a coordinate which runs along the panel. The second 
point is associated with the applicability of the in-plane boundary conditions at the 
element level and hence the evaluation of the nonlinear stretching forces. Gener- 
ally speaking, the nonlinear stretching force is constant along the panel when the 
in-plane inertia forces are neglected, and there are no distributed in-plane loads. 
For a two-dimensional panel with immovable edges this force at the system level is 

ra 

proportional to the integral / (w, x f dx , which is evaluated along the panel length. 
However, at the element level this is no longer true, since a term related to u, x must 
be considered. Neglecting to recognize these in-plane displacementes leads to the 
assumption that they are zero at all nodal points, and that the nonlineai stretching 
force varies along the panel. 

In the linear case [65l-[73] the analysis is very much similar to the one used for 
Galerkin’s method, equation (1), with the difference that the matrices [M], [ I \ ], and 
[A] for the system are arrived at by assembling individual matrices for each element. 
Also, the generalized coordinates {9} represent now the nodal displacement (s). The 
buckling phenomenon due to in-plane stresses arising from initial and/or thermal 
effects can be taken into account by introducing the concept of geometrical stiff- 
ness 174 '. In the nonlinear case [75 ' _l84 ' the eigenvalues for a given dynamic pressure 
specified by A must be determined iteratively. The most used iterative procedure is 
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explained in detail in the work by Mei and Rogers Briefly speaking, for a given 
A the linear flutter problem is solved first, yielding a first approximation for the 
mode shape, which is normalized by its maximum component. Then the first 

approximate value for the displacement is given by {q}i = c 5R( {<^} 0 e iu; * ), where 
c is a given amplitude of panel oscillations, and 5R( ) denotes the real part. It is 
thereupon possible to evaluate the nonlinear stiffness matrix, which depends on the 
displacement, and to use this matrix to obtain a new approximation to the mode 
shape, associated with the amplitude c. The iterative process is repeated until 

a convergence criterion is satisfied. In many cases it is most efficient and accurate 
to base the convergence criterion merely on displacement quantities, such as the 
norms defined by Bergan and Clough 


1. Modified absolute norm: 
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Here N is the total number of unknowns components, and Ar, is the change in the 
displacement component i during a given iteration cycle. Every such component is 
scaled by a reference displacement quantity r, ire f, which are, in general, not equal to 
the corresponding total components because if r, is close to zero, the ratio Ar^/r, 
could be a large number even after convergence has occurred. Instead, every Ar, is 
scaled by the largest displacement component of the corresponding “type” (deflec- 
tions/rotations). A frequency norm, defined as the absolute value of the ratio of the 
change in eigenvalue during a given iteration cycle and the eigenvalue itself, may 
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also be useful in the convergence decision. Once the converged values for the eigen- 
values are known, the problem is solved as in the linear case, that is, \ CT corresponds 
to the lowest value of A for which coalescence occurs (usually for c = 0). 

Another possible method of analysis is the so-called second method of Lia- 
punov 1^1. It has proved useful in a number of stability problems involving ordinary 
differential equations, but the first application to a panel flutter problem involving a 
partial differential equation is believed to be given in the work by Parks Briefly 
speaking, the Liapunov method involves finding a functional V (defined in phase 
space) which must possess certain sign properties as prescribed by one of a number 
of stability and instability theorems an( j that describes the motion of the panel 
when it is disturbed from its equilibrium position. For example, if V is a positive 
definite function and dV/dt can be shown to be a negative semidefinite function 
along every trajectory, then V — ♦ 0 as t — > oo and the disturbed motion must die 
out. The fact that such a function cannot be found does not imply that the system 
is unstable. Indeed, the main drawback of the method is that there is no established 
procedure for producing a Liapunov function for any given dynamical system. 


1.3 Present Application 

The present application addresses the issue of the breakdown of continuum fluid 
mechanics when considering panel flutter in the case of hypersonic flight vehicles 
traveling at high altitudes, as would occur along the upper trajectory of the National 
Aerospace Plane. Once the atmosphere becomes very rarefied, one can expect such 
a breakdown, especially for surfaces parallel to the free stream or in the leeward 
part of the vehicle. The simplest possible way to deal with the aerodynamics in 
this case is to neglect the effect of collisions among molecules, which corresponds 
to assuming that a free-molecule flow exists. Another approach is to use some kind 
of “particle simulation method” to take the effect of molecules intercollisions into 
account. The most important drawback in the latter case is related to CPU time 
requirements, even on the largest supercomputers, as well as memory limitations. 
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As of today, such applications do not seem to be feasible, although they are possible 
from a conceptual point of view. 

In this dissertation the free-molecule airloads are used in the study of panel 
flutter in the case of a two-dimensional (beam-like), simply-supported panel. The 
equations of motion are derived in detail in Chapter 2, allowing for the existence 
of distributed longitudinal loads, as well as distributed bending moments, along 
the panel. These kinds of loads are necessary in order to model the aerodynamic 
shear effects. Linear springs in the longitudinal direction at the panel’s leading and 
trailing edges are also introduced in the structural model. Neglecting rotary and 
longitudinal inertias, one can reduce the problem to an integro-differential equation 
in terms of the transverse displacement. The problem is solved using Galerkin’s 
method and direct numerical integration in the time domain. 

Chapter 3 deals with the aerodynamic models. The first part consists of a 
brief summary of existing formulations in the case of continuum flow. Then the 
free-molecule airloads are obtained using a quasi-steady approximation. This ap- 
proximation consists of assuming that the angle of attack at a given point in the 
panel is equal to the local slope, w, x , plus the angle induced by the panel nor- 
mal velocity, w, t /U^. It turns out that the aerodynamic shear stress, considering 
the hypersonic regime, is at least one order of magnitude larger than the pressure. 
Finally, particle simulation methods are briefly discussed. 

Chapter 4 contains the presentation and discussion of numerical results, with 
the response of the system being analyzed for a wide range of parameters. Finally, 
Chapter 5 is a summary of conclusions and discussions of possible future work. 
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STRUCTURAL MODEL 


2.1 Equations of Motion 

As discussed in the Introduction, the only case to be considered in this work cor- 
responds to a two-dimensional panel (beam-like panel) with some small initial cur- 
\atuie. A consistent way to obtain the equations of motion in this case is provided 
bj Steele^, the effect of moderate rotation of the panel midplane surface being ac- 
counted for. Actually, the static version of the governing equations in the case of a 
flat panel are those of the von Karman plate theory. It should be pointed out that 
Kirchhoff’s hypothesis is assumed to be valid, implying that the effect of the trans- 
verse shear deformation is neglected. Roughly speaking, the length/thickness ratio 
of a isotropic panel should not be less than 15 if results with reasonable accuracy 
are to be obtained [91 k In the case of composite panels the previous number should 
read 25. 

A complete derivation of the equations of motion is given in this section. Con- 
sider then a panel of length a as shown in Figure 1. The displacements from the 
unstressed state of the panel’s midplane surface in the x and 2 directions are de- 
noted by u and w, and the total transverse displacement of a given midplane surface 

,This reference corresponds to part of the classnotes provided by Prof. Charles Steele for the 
course ME 241 at Stanford University, CA (Moderate Rotation Theory for Beams and Shallow 
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z, w 

Pz 



Figure 1: Two-dimensional panel with initial curvature, 
point after deformation is given by 

^totaK^M) = w 0 (x) + w(x,t) ■ (4) 

In equation (4) te 0 (ar) indicates the initial undeformed shape of the midplane sur- 
face, while w(x) corresponds to the transverse displacement of the midplane surface 
relative to its undeformed configuration. 

The strain e of the midplane surface in the x direction is given by 

£ 'Uix + Woix iv ,x d* ^ n? , x . (5) 

Note that the subscript ( ), x denotes differentiation with respect to x. Also, the 
terms in w in the equation represent the additional stretching of the midplane 
surface due to moderate rotation. 

When following the procedure laid down by Steele (see footnote in page 12), the 
internal forces in the panel’s cross-section are defined as the axial stress resultant 
A^, the bending moment A/, and the shear force Q. These resultants, with sign 
conventions defined as in Figure 1, can be expressed as 


N 

= Eh 

[e-aAT(x)] = 



= Eh 

[ (iz, x + u> 0 , x w, x + \ w 2 , x ) - q AT(.r )] , 

(G) 

M 

= Dk 

= - D u>, X£ , 

(7) 

Q 

= M, x 


(S) 


Here a is the coefficient of thermal expansion, A T(x) is the temperature difference 


CHAPTER 2. STRUCTURAL MODEL 


between the panel and the supports, h is the panel thickness, D is the panel stiffness, 
defined as 

D = 12(1 - v 2 ) ' 

and k is the curvature change of the midplane surface, defined by 

K = XV , xx 

All other symbols are defined in the Nomenclature table. 

The distributed loads p z and p x being applied to the panel can be represented 


Pz = Pz ~ Pz 


Px = Pt-Px » ( 10 ) 

where the superscripts ( ) + and ( )“ correspond, respectively, to the upper and lower 
surfaces of the panel. It should be noted that the sign convention for pf and p* 
is the usual convention for stresses. Then pj corresponds to the stress a zz on the 
lower surface and its positive direction is toward the negative 2 direction. With the 
above representation for p x , when p+ / -p~ there must exist a distributed bending 
moment, m x , along the midplane surface. If rotary and longitudinal inertial loads 
are neglected, and the panel is thin, it is reasonable to take 

m x = - {pt + Px) ■ ( 1 1 ) 

Finally, the distributed loads, with the sign conventions illustrated in Figure 1, can 
be expressed as 

P , = + + . < 12 ) 

Pi = Px(Zyt) + APi tal ( x ) > 

m x = | [pJ(x,() + Apy(*)] . (14) 

The first term in equation (12) corresponds to the transverse inertial load, while 
the superscrit ( ) A indicates an unsteady aerodynamic load and the superscrit ( ) 
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indicates some kind of static load. The kind of aerodynamic load that gives rise to 
P x a shear stress. It is assumed that there is no shear stress in the lower surface of 
the panel in order to write m x as in equation (14). This last assumption corresponds 
to the panel being exposed to airflow in its upper surface, and the cavity behind it 
being filled with fluid at rest. 

The equations of motion, as well as the appropriate boundary conditions, can 
now be obtained from Hamilton’s principle: 

/ (N $£ + MSk) dx — f ( p x 6u — p z 6w + m x 6w, x ) dx + 

Jo Jo 

+ ki u(0, i) £u(0, t) + k 2 u(a, t) 6n(a, t) + 


- (N m 6u + Q*6w-M*6w, x ) 


a 


= 0 . 
o 


(15) 


Here N m , Q *, and M* are the end loads acting at x = 0 or x = a, while k, and k 2 are 
the spring constants of longitudinal springs attached at the panel ends [the subscript 
(•)i corresponds to x = 0, with (•) 2 corresponding to x = a]. It should be emphasized 
that the first integral term in (15) corresponds to the variation of the strain energy 
of the panel, while the second integral term is associated with the work done by 
the external distributed loads. The potential energy of the longitudinal springs is 
taken into account by the terms in the second line of (15). Finally, the work done 
by the end loads is associated with the terms in the third line of equation (15). 

Using the appropriate expressions for e and k, and integrating by parts, one can 
write the variation terms in the integral portions of equation (15) only in terms of 
6u and Sw. These steps give rise to the following equations of motion: 


N, x + p x = 0 , 


(1C) 


M, xx + [N(w 0 + w), x ]^+ p z + m T , x =0 . (17) 

Equations (16) and (17) can be seen to correspond to equations (1-54) and (1-5G) of 
Reference [91], pp. 17, if these are conveniently reduced for the case of cylindrical 
bending of a flat panel. It is worth mentioning that the equations in the latter 
case were derived from the nonlinear theory of three-dimensional elasticity, with 
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the basic assumptions that: (1) linear strains and the squares of angles of rotation 
are small compared to unity; (2) the rotation about the normal axis to the panel can 
be neglected in the calculation of strains and stresses; (3) Kirchhoff’s hypothesis is 
valid; and (4) the material is linearly elastic. The appropriate boundary conditions 
at x = 0 or x — a are 


N = N* ± k{U or u = u* 


M = M* or w, x = w‘ x 


M, x + N(w 0 + w), x + m x = Q* or w = w* 


(18) 

(19) 

( 20 ) 


The upper sign in front of the ku term in the forced boundary condition (18) 
corresponds to x = 0, while the lower sign to x — a. 

From equations (16) and (17) it is clearly seen that the resultant force N must 
be written in terms of w in order to render the equation of motion in the transverse 
direction complete in itself. By assuming that the distributed load p x is known, 
equation (16) can be integrated with respect to x from a generic value x to x = a, 
yielding 

N(x,t) = N(a,t) + [ p x (T),t)dr) . (21) 

J X 

Also, recall that at the boundaries 


N = N m ± ki u 


(18b) 


Then, if x is taken equal to zero in equation (21), one can obtain a relation between 
u(0, t) and it(a, t ): 


u(0,t) = 


*i 


N*\ x=a 


N* 


[x=0 


+ / Px(x, t) dx — k 2 u(a, t) . (22) 

Jo 


The next step is to integrate equation (6) from x — 0 to x = a. It is going to be 
assumed that the panel is uniform, and u>(0,<) = w(a,t ) = 0. Using the previous 
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results for N(x,t ) and u(0,t) gives 

Eh r a 


» ▼ / ,, Eh f a 2 Eh f a 

N(a,t) = a k — / w, x dx~a k / ww 0 , xx dx + 

Za Jo a Jo 

f a f a 

/ / p x (r],t)dr} dx 

Cl Jo J x 


+ <* k 


a ~ I p x (x,t ) 

Or*:, Jo 


+ 


dx + N x . 


Here 


= 


&ko = 


a k = 


k x a 


k i ci “f* E h 

k 2 a 

k 2 a -f- Eh 

^A:2 

a /:i + a k 2 ~ a k x <*Jfc 2 


(23) 

(24) 

(25) 

(26) 


N, = (1 - — ) N'\„ a - a t N ' 

Eh 


x = 0 


+ 


Eh r a 

— ot k a / A T(x)dx. (27) 

a Jo 

Note that u(0,t) = 0 is achieved if a k} = 1, that is, k l = oo. In this case a k = a kl . 
Similarly, a ki = 1 leads to u(a,t ) = 0 and ot k = <x kj . From the definition of N x , 
which corresponds to the net applied force at the panel’s ends, one can clearly see 
the influence of different temperature distributions along the panel, A T(x). Given 
two temperature distributions A T x {x) and A T 2 (x), the effect on the panel response 
is identical if 

r ATi(x) dx = r AT 7 (x)dx . 

Jo Jo 

This fact has been previously mentioned in the work by Xue and Mei^l], where 
the finite element method was used. 

The equation of motion in the transverse direction is found to read 


Dw, xxxx + p x (u>o + rc), x — |-/V(a,f)-f J p x (rj,t)di ] 

+ Pph w, tt - p * - m x , x = A p 


(u ’ 0 + w ) , xx + 


stat 


(28) 
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It is easily seen that, when p x = 0, as well as u> 0 (x) = 0, equation (28) reduces to 
the form used in classical panel flutter 149 '. In this case it is worth emphasizing that 
the resultant force N is constant along the panel at any given time, contrary to the 
present case. 

One may now assume that both Ap* tat and A p* tat are constant along the panel. 
If only a flat panel is to be considered, i.e., w 0 (x) = 0, equation (28) becomes 

Dw, xxxx + p p h w, tt + [pxOM) + Apf 1 ] w, x -p x - ~(p x ),t + 



Ok 

a kl + 2 





+ 


Ok_ 

a 


f / P x (r],t)dr]dx + a k — / p x (x,t)d. 

JO Jx &k\ J® 


X + 


+ 


I Px(iht)drj 

J X 



W , 


= Ap 


s tat 


(29) 


The variables in equation (29) are made nondimensional as follows: 


i = - ■ w = j , 

a n 


(30) 


T = t 


D 


Ppha 4 


The result reads 


W,a'' + W, TT + P*(t,T) + P x W*-pf - f 


6o fc (1 - v 2 ) W% d v + (l - ^ ^ - *) P x + + 

- a k [ f Pz(r],T)dr]dZ + a k — f p x {r},T)dr] + 

Jq J £ Qk\ 

+ J Px{p,T) d v + Ri 


= P' • 


(31) 
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The following definitions were also used in equation (31): 


P A X 

II 

H^> 

P x = 

P r, 

Pz 

qA a 

= hD P -' 

Pz = 

a 4 

A »* tat 

hD lz 


R x 

II 



(32) 


2.2 Initial Curvature or Imperfections 


A quick look at equation (28) may suggest that only the explicit terms in w 0 have to 
be included in equation (29) in order to obtain the equation of motion which takes 
into account the effects of initial curvature. However, a more careful examination 


shows that the inherent geometric curvature of the panel also induces steady 
dynamic loads, which should be included into A ;/ tat (*). If the aerodynamic 


<UTO- 

l()rl<ls 


are assumed to be of the form 


P*(Z,T) = AiW,z+A 2 W,t , ( 33 ) 

Jj£(£,T) = A 3 W < +A 4 W, t , (34) 


one can write 

Ap* tal (0 = P x + A 3 W 0 , € , ( 33 ) 

A pf at (0 = Pz + AWo* • ( 3 °) 

Note that the present definitions of P x and P z are similar to the earlier ones, le- 
taining only the constant part of the Ap stat load along the panel. 
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At least two kinds of initial curvature may be of interest: 


Constant curvature panel : such a panel was studied by Dowell t 52 - 53 ) j n the 

late 60’s. The initial geometry of the panel may be approximated by the 
parabola 



H 

h 



(37) 


where H is Ihe maximum height of the curved panel. Then the nondimensional 
initial curvature of the panel is just 


r 0 = -W 0 , u = 8 ^ . 

h 


(38) 


2- Panel with sinusoidal curvature : this case corresponds to take 

wo H . 

Wo = ~h = ~h smpn t ’ (39) 

where H is again the maximum height of the curved panel. Then 

H 

Po = -^o,« = ~r (p*) 2 sinp-n-Z . (40) 

Such a curvature may be of interest when approximating a general initial 
geometry using a Fourier series. 


2.3 Role of Structural Damping 

Since this work deals with a rarefied atmosphere, the magnitude of the aerodynamic 
damping terms, with coefficients A 2 and A 4 , may be expected to be very small. This 
raises the question of how to include more damping, if necessary, into the equation 
of motion. During numerical integration such damping may be required to yield a 
converged solution. 

Structural damping can significantly modify the panel flutter boundaries, but 
the modification is extremely dependent on the type of structural damping model 
employed. If only linear damping is considered, the work by Ellen [92) provides a 
useful classification of this mechanisms and its influence on the flutter boundaries. 
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Basically, structural damping can be introduced into the system by adding a term 
of the form 

d }+1 w 
^ dt dx 3 

to the equation of motion (28), where g is a structural damping coefficient. Two 
possible mechanisms can be classified as follows: 

1. Viscous damping: in this case g is a constant. In its simplest form, j = 0, 
the term is of the same type as the aerodynamic damping term in the case 
of linear piston theory, as well as the present rarefied formulation. This has 
the analytical advantage of grouping the structural and aerodynamic damping 
terms in a single total damping quantity. This type of damping is always 
stabilizing. However, when j > 0 the system can be destabilized since the 
structural damping term may then supply energy to the system instead of 
dissipating it. 

If a modal solution such as Galerkin’s method is used, damping may be taken 
to vary from mode to mode; this procedure, originally with j = 0 terms, is 
equivalent, on varying g , to the insertion of an appropriate j > 0 term with 
constant g. For example, in Reference [36] a two- mode solution is used to 
investigate the effects of structural damping with C1/C2 = 1 and 4. Here £ n 
is the damping ratio of the nth mode and g = 2£ n <^o being a suitable 

reference frequency and u n the nth natural frequency. Note that C1/C2 = 4 
implies Ci ^1 = ( 2 ^ 2 , that is, g = constant whence j = 0. In this case the 
critical flutter parameter increases with structural damping, although this 
effect is almost not noticeable for small values of aerodynamic damping, <j,\ 
On the other hand, the case C1/C2 = 1 can be shown to correspond to j — 2, 
and the effect of increasing the structural damping is destabilizing, especially 
when g \ is small. Lottati studied the case for j = 4 assuming a viscoelastic 
type of damping, which corresponds to multipling the panel stiffness D by a 
factor equal to 1 + g — . His conclusion is that this kind of damping has a 
strong destabilizing effect for all values of g ± 0. One point that is worth 
emphasizing after all this discussion is that the assumption of = Ci^i? 
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Figure 2: Correlation of predicted damping ratio for j = 0 with Zener equation. 

corresponding to j = 0, yields values of C„/Ci that follow closely the results 
obtained if the “Zener equation” described in Reference [94], 

U / CJ re l a x 

1 + (u;/u> re lax) 2 

is used for values of frequency at least four times larger than the relaxation 
frequency u> re i ax defined there, see Figure 2. This may be significant since 
the Zener equation is based on a very basic physical assumption, i.e., that the 
mechanism which causes energy dissipation in metals is heat flow due to strain- 
induced temperature gradients. Experimental data obtained for aluminum 
panels correlate well with the theoretical results. 

2. Hysteretic damping: it can be shown that the viscous type of structural damp- 
ing just discussed results in a dissipation of energy per cycle which is piopor 
tional to the frequency of oscillation. As structural hysteresis loops appear 
experimentally to be frequency independent, the previous damping expres- 
sions need to be modified to include this mechanism. This can be done by 
taking g = < 7 /|u>|, where g is a constant parameter and |u| is the modulus of 
the complex frequency of the motion, which makes this definition applicable 
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to nonharmonic oscillations. The drawback of this kind of formulation is that, 
for a real frequency of vibration, this term becomes i g d^u> j dx^ for the positive 
frequency spectrum, a form that is strictly valid only at neutral stability. 


2.4 Galerkin’s Solution 

As discussed in the Introduction, one of several ways to solve the mtegro-differential 
equation (31) is to use Galerkin’s method, where a expansion of the form 

OO 

W(Z,T)= J2 <lm(T)W m (0 (41) 

m~ 1 

is adopted. Here W m (£) are assumed modes and only a few modes are used in 
practice, that is, the series is truncated. In the case of a simply supported panel it 
is possible to take 

Wm(0 = sinm7i-£ . (42) 


By applying Galerkin’s method, one converts the integro- differential equation (31) 
into a system of ordinary differential equations in time. The number of equations 
equals the number of assumed modes taken in equation (41). For the case of modes 
(42), the nth equation is given by 

a n 9n + ( a nr Qr + « nr <JV + ^*** <7r) Qn + 


X/ ^ nm ^ ^ ^ ^ nmr 4r T ) ] C nmr Qr J Qm T 

m m \ r r J 


^ v ^nm Qm Qn ~ I - ^ ^ ^ nm Qm 4“ ^ Qn — * (43) 

m m — ' 

The definitions of the various coefficients are listed in Appendix A. Note that the 
introduction of viscous structural damping with j = 0 in the present case amounts 
to add the term 


(44) 


to the term c„ q n . 


K 2 Cl <ln 



Chapter 3 

AERODYNAMIC MODELS 


3.1 Continuum Flow 

In the study of panel flutter in this regime, the aerodynamic shear loading is usually 
neglected compared to the aerodynamic pressure. Also, the latter may be considered 
as the sum of two parts: 1) one given by the pressure fluctuations on the panel in 
the absence of any panel motion, e.g., due to turbulent boundary layer fluctuations, 
and 2) the other due to the panel motion itself. Inherent to this superposition is the 
assumption that the panel motion and the consequent portion of the aerodynamic 
pressure are sufficiently small so that the turbulent pressure fluctuations themselves 
are not substantially modified This means that this last contribution can be 
taken as prescribed from experiment or separate analysis. On the other hand, the 
aerodynamic pressure due to the panel motion may be related to w in general by 
an expression of the form (two-dimensional case): ^ 

+ f [ A(x — x',t — t') w , x < + jt— w ,(> dx'dt'\ . (45) 

Jo Jo L Ooo J ) 

Here A is an “indicial” function which depends parametrically on the freestream 

Mach number, = Uoo/a.oo , where a <*> is the freestream speed of sound. Note that 

this term describes the effects of spatial and temporal memory, i.e., the pressure at 

24 



3.1. CONTINUUM FLOW 


25 


a particular point r, at a particular time is influenced by the motion at all points, 
0 < x' < a, at all previous times, 0 < t f < t. The neglect of the memory effect leads 
to the so-called linear “piston theory” approximation $5»4] 

p * (x ’ t) = ~ e i£r h' + [£:”•'] = (4&) 

Poo^oo \U 00 w^ x -\- w ] . (46b) 


This result is accurate enough for large values of M ^ (high supersonic-hypersonic 
regime) and can be recognized as the relation between the motion of and pressure 
on a piston in a tube. The total piston velocity includes both a convection term, 
Coo as well as the direct velocity, xv,(. Here the panel is the equivalent piston 
and the fluid “tube” is perpendicular to the undeformed panel. 

By nondimensionalizing according to expressions (30), one modifies equation 
(46a) into 


p?({,T) = - A 


W* + 



(47) 


where: 


2ga 3 _ Poo a 

Moo D ^ p p h 

v = \p~ui . 


(48) 


Note that the parameter p can be identified as a mass ratio. Also, q is the familiar 
dynamic pressure. 

Another level of approximation can be obtained if simple harmonic motion is 
assumed and equation (45) is expanded as a power series in frequency, u>, that is, 


pf(x,u) = - 


P~Ul 


\fMi 


+ 


ML - 2 


ICO 


M 2 

oo 


1 u* 


W 


+ 0(uj 2 ) 


(49) 


In the time domain this corresponds to the so-called “quasi-steady” approximation, 
and it is clear that it reduces to the linear piston theory result for 1. At 

low supersonic Mach number the second term between the brackets gives rise to 
negative damping when < 2, which can lead to single-degree-of-freedom flutter. 
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All the remarks regarding the applicability of piston theory results also apply to 
this approximation. 

For the sake of completness, it should be noted that second-order piston theory 
introduces nonlinearity in the aerodynamic model by adding the term 



to the linear expression (46a). This nonlinearity may become important for hy- 
personic speeds because of the Mach number factor multiplying the term between 
parentheses, but Eastep and McIntosh emphasize that only the terms 

_ ^ M °» [(•»» )’ + 7T 1 } (50) 

Moo 14 l Voo J > 

are relevant for parameter ranges anywhere near those found in practice. In fact, 
the term proportional to (w, x ) 2 produces an overpressure, tending to push the panel 
into the cavity, for any excursion of the panel from its flat undisturbed state. 

If “external” aerodynamic pressure contributions independent of the panel mo- 
tion are not considered, and linear piston theory is used, the values foi the -4, 
coefficients in equations (33) and (34) are 



a 3 = o, A 4 = o. 


3.2 Rarefied Flow 

As a spacecraft or aircraft gains altitude, the air becomes increasingly less dense. At 
a sufficiently high altitude the atmosphere becomes so rarefied that one can expect 
the breakdown of the continuum hypothesis. When this happens, the discrete molec- 
ular nature of the gas surrounding the spacecraft can no longer be ignored. At this 
point methods from kinetic theory t 96 H 98 l must be used to predict the aerodynamic 
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behavior. It is useful to mention that the analysis of such flow fields is included in 
the field of rarefied gas dynamics. The available literature in this subject is very 
extensive, but the proceedings of the biannual series of International Symposia on 
Rarefied Gas Dynamics form a unique record of the field. More especifically, the 
recent survey by Muntz is a very enlightening work. 

The basic parameter for specifying the degree of rarefaction of a low-density 
flow is the Knudsen number, Kn, which is the ratio of the molecular mean free path 
of the gas to a characteristic dimension of the flowfield. Continuum flow theory 
is assumed to apply when this parameter is much smaller than one. On the other 
hand, a large Knudsen number (usually Kn > 1) characterizes the so-called free 
molecule flow in which intermolecular collisions are neglected. The region between 
these limits is generally referred to as the transition flow regime. 

3.2.1 Free-Molecule Regime: Quasi-Steady Approximation 

The estimate of the aerodynamic loads due to the panel’s motion in the free- molecule 
flow regime is based on the previous development made by Ashley In this study 
the steady values of the pressure and shear stress exerted on a surface element at 
an arbitrary orientation with respect to the mean flow were obtained through a 
general procedure. Note that only a simple gas (i.e., the gas is composed of a single 
chemical species) is considered. Assuming | w, x | <C 1, one sees that at any instant 
the “steady” slope of a surface element of the panel to the mean stream 17, X( is 
w, x . The “quasi-steady” approximation implies that the surface element normal 
velocity w, t induces a further inclination relative to the mean flow, such that the 
total slope of a surface element becomes w, x + w, t / U oo- Inherent to the use of this 
approximation in the procedure described by Ashley is the assumption that enough 
molecules strike the surface element during an interval At, and that both w, x and 
w, t do not change appreciably during this time interval. 

Following the development by Ashley, the normal and shear stresses are sepa- 
rated into those due to the incoming stream of molecules (p,- and r t ) and those due 
to the reflected or reemitted stream (p r and r r ). The total pressure and shear stress 



28 


CHAPTER 3. AERODYNAMIC MODELS 


are then given by 


T = Ti + T r . 
P = Pi + Pr , 


(52) 

(53) 


The incident values can be obtained by considering first a reference frame x'y'z 1 
fixed relative to the mean motion of the gas. Then the number of molecules of class 
c' (i.e., molecules having velocities in the range c' to c' + dc ') per unit volume is 
given by ^8] 


dn — Hoq/max dc\ dc' 2 dc 3 = 

3/2 


= Tlr 


//?oo\ e -(}oo(c?+c?+c‘ 0 dtc\dc' 2 dc' 3 , 

\ K / 


(54) 


where = Poo/ rn is the number density of the flow, is defined by 

1 


/? oo = 


2RT n 


m is the mass of one molecule, and R is the gas constant. Note that the Maxwellian 
distribution of thermal velocity, /max, was used to characterize the velocity distri- 
bution function of the undisturbed flow. 

Consider now a surface element moving with velocity components e i , fi 2^001 
e^Uoo in the directions x', y', z' respectively. A new reference frame x a y s z a is instan- 
taneously fixed to the surface element, parallel to the x'y'z' system and such that 
x a is the outward normal to the surface, as shown in Figure 3. Then any molecule 
has velocity components in the x s , y s , z a directions given by 


c' -eU 0 


(55) 


Hence the distribution of molecular velocities, which was initially given by equa- 
tion (54), can be rewritten in terms of Cj, c 2 , and c 3 : 

3/2 


dn = rioo (— ) e -/ 3 °°[( c ' +e i t/o °) 2 +( C2 + e2l/oo ) 2+ ( c 3 +e 3 l/o o) 2 ] dc x dc 2 dc 3 . 


(56) 




Figure 3: Reference frames for the quasi-steady approximation. 


The flux of molecules of class c arriving at the surface element is given by 
— Ci dn ^1, that is, 

™oo ci/max^ci dc 2 dc 3 . (57) 

The minus sign in the above expression is necessary because only molecules with 
velocities in the x s direction in the negative range — oo < Ci < 0 can strike the 
surface. If this expression is integrated in each velocity range, one can find the 
number of molecules striking unit area of the surface element per unit time: 


Mi 


/a \ 3 / 2 o 

™oo f — ■ ) f Ci e“^ 0 °( Cl+ei ^ 00 ) 2 dc x x 

y 7T y J — oo 


/ oo „ roo 

e -^(c 2 + e 2 t/oo) dc 2 r e -3.Mc, + c,V.^)- (jci _ 
- oo J - X ' 


n r 


e I U 




+ 


1 + erf 


l ^ oo yfc) 


Jl JJ 

,l OO OO * / \ 

^a(ei) • 
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(58) 
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In equation (58) the following functions were used: 


* 3 (0 = 

(0 H — $ 2(0 , 

s 

(59a) 

*i(0 = 

-c 2 

e 

y/Zs ’ 

(59b) 

<M<0 = 

1 -(- erf (c) . 

(59c) 


Here erf(c) represents the conventional error function of statistical theory, defined 

by 

2 2 

erf (c) = — — / e n drj . 
v/?r Jo 

The parameter s = Uoo/c mp is the so-called molecular speed ratio or the “Mach 
number” based on the most probable molecular speed, c TOp = 1 j \J . It should be 
noted that for a perfect gas s = \Jyj2Moo. 

The forces on the surface can be found by considering first the total momen- 
tum / of the mass stream mA r , in some arbitrary direction specified by the di- 
rection cosines £ t , £ 2 , £3. This is found by integrating the momenta of molecules, 
— ( c i £\ + c 2 £2 + c 3 ^ 3 ) c x dn, in each velocity range, yielding 

= — <1 |(<h £\ + e 2 £2 + C3 £3) $1 (ej ) + 

+ + c l ( e l £\ + e 2 4- f'3 £ 3 ) ^ 2 (^ 1 ) | • (GO) 

The pressure p t is then obtained by setting £ x = — 1, £ 2 = £3 = 0, while the shear 
stress r, corresponds to £ 2 = 1, £\ = £3 = 0. 

In the case of a two-dimensional, rigid surface inclined with respect to the undis- 
turbed flow by an angle #, one has c x = sin 8, e 2 = — cost), and e 3 = 0. For the cor- 
responding panel flutter case it is reasonable to take \8\ <C 1, such that sin# 8 and 
cos 8 rs 1. The quasi-steady approximation then calls for setting 8 — w, x + w, t /U^, 
where it is important to point out that w, t /Uoo < 1 because the flow is considered 
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to be hypersonic. With these results, the incident pressure and shear stress are 
given by 


r, = g$ 3 (e) , 

Pi = Q {^$3(0 + 2^2 $2 ( c )} ’ 


(Cl) 

(62) 


where 


e = s 


(-♦a) 


Note that in the case of a steady tangent stream (w, x = u>,< 
equations reduce to 

Q Poo U(X> 


‘steady 0F s 2 


■steady 2 S 2 4/? 0 


= 0) the previous 

(03) 

(04) 


Before turning to the problem of finding the reemitted pressure and shear stress, 
it is important to comment on the gas-surface boundary conditions. Despite a great 
deal of theoretical and experimental study on this subject U00.103]^ there is still no 
general consensus about what kind of model should be used when tieatiug gas- 
surface interactions. The two most simple models are those of specular and diffuse 
reflection ^4) 

The specular model assumes perfectly elastic collisions between molecules and 
the surface, that is, the molecular velocity component normal to the surface is 
reversed, while that parallel to the surface remains unchanged. Thus 


Pt spec 


P. • 


( 60 ) 


This is a useful model for analytical studies but it is not physically meaningful in 


real gas-surface interactions. 
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In diffuse reflection the velocity of each molecule after reflection is independent 
of its incident velocity. In fact, it is assumed that the molecules are brought to 
rest relative to the surface and then reemitted with the equilibrium distribution 
corresponding to a temperature T r , which may differ from the temperature T p of 
the surface. Thus the values of T r and p T correspond to expressions (63) and (64) 
applied for a stationary gas at temperature T r , that is, 


r 7diff = 0 ’ 

p r n T m 

Prdi * = = 4&~ ' 


(67) 

(68) 


In the absence of adsorption or emission effects at the surface, the incident number 
flux of molecules A", to a surface element must be balanced by the reflected number 
flux Af r . Therefore, using equation (58), 


n r 




2 


* 3 (e) 


n T 

2 y/nJTr 


=> P.dilf = TTJ, $3(e) . 

Actually, the factor multiplying the function $ 3 in the above expression can be 
reduced to a more suitable form if the definition of the speed of sound for a perfect 
gas is used, and the relationship between and s is recalled. Doing this 


P 


r diff 



(69) 


Note that in the case of a steady tangent stream the expression for becomes 


= i [TL - p °° [TL 

Pr ' steady 7 Ml V ~ A(3oo V Too 


(70) 


The extent to which the reflected molecules have their temperature adjusted 
toward that of the surface may be indicated by the “energy accommodation coeffi- 
cient” a e , defined by 

Si - Sr 


( 71 ) 
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In expression (71) £, and £ r are respectively the incident and reflected total energy 
fluxes, while £ p is the total energy flux for the case of complete diffuse reflection, 
i.e., T t = T p . Some experiments with engineering surfaces in contact with gases indi- 
cate that the reflection process approximates complete diffuse reflection. However, 
Hurlbut mentions that recent measurements made in the Shuttle Orbiter, as 
well as the study of the decay of rotational frequency of the spin stabilized satellite 
Explorer VI, rule out the complete diffusive process as the predominant scatter- 
ing mode of molecules interacting with surfaces. Alternative empirical reemission 
models have been proposed in past years. An example is the Nocilla model studied 
by Hurlbut and Sherman [10o] , as we jj as Pandolfi and Zavatarro^^G]. The inclu- 
sion of some arbitrary fraction of specular reflection appears, however, to remain 
the most practical way of allowing for departures from complete diffuse reflection. 
This approach is going to be followed here, such that a fraction c\ ni of the molecules 
is considered to reflect specularly, while the remaining fraction (1 — n m ) reflects 
diffusively with T r = T p . Note that the coefficient o m can be viewed as a “momen- 
tum accommodation coefficient”. According to the foregoing convention, and using 
equations (52) and (53), the total shear stress and pressure exerted on unit area of 
the surface element can be written as 


r = ( 1 - a m ) r, , 


P — (IT ) Pi T ( 1 0 7)l ) p r< 


mr 


With the use of equations (62), (61), and (69) the final expressions for r and p are 

- = (l-a m )$ 3 (0 , (72) 

<7 


P 

q 


/ -i . \ C 1 a m / ^ 7 ), 

( + ° m) s + V2 7 


*1* a ( ^ ) H — T'if ( ) • (73) 

Is* 


These relations are still nonlinear in the motion but they may be lin- 

earized if one recognizes that the parameter e may be reasonably taken to be small. 
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Expanding the exponential and “erf” functions in and $ 2 , assuming that e < 1, 
and retaining only terms of up to first order in e gives 


r 

9 

P 

9 


1 ^ m 

\Ztts 


+ (1 — «m) 


U>, x + 


— 

u a 


1 

2 s 2 


1 + + (1 
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(74) 


+ 



4(1 + a m ) + 7r(l - a m ) 




+ 


17 . 


(75) 


Note that both expressions have a steady part that does not depend on the panel mo- 
tion, as should be expected. Also, the ratio of p/r for both the steady and unsteady 
parts is of order of 1/s. Since the regime considered is hypersonic {M^ > 10) and 
M^/s = 0(1), the results indicate that the pressure is much smaller than the shear 
stress. It must be pointed out that equations (74) and (75) correspond to equations 
(2.55) and (2.56) of Reference [107], pp. 43, where the parameter 9=1 — a m . 

In order to find and ]/} it is necessary to decompose both r and p along the 
x and axes of the structural model, and take only the unsteady part. Thus, 



Note that the term [(• ••)/s 2 ]iy, x due to the contribution of the steady part of the 
pressure p to the load p £ has been neglected. On the other hand, the corresponding 
contribution of the steady part of the shear stress r to p * was retained because it 
is of the same order of magnitude as the other terms due to the unsteady part of p. 
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If equations (76) and (77) are nondimensionalized according to definitions (30) and 
(32), the final result reads 




A M n 


th (1 - « m ) 


W < + 


A M 0 


W, 
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(78) 
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2(1 + a m ) + £(l -a m )Ve 
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W T 

. 
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(79) 


where: 

th = - , e = 2k. . (so) 

a 1 oc 

Note that three more parameters are added to the problem when t h(' effects ot shear 
stress and rarefied flow are introduced: th, a m , and (~). In fact, th is inherently built 
into the definition of P x if it is to be compared to P z , and the ratio T v /T (yj can be 
taken as a variable of the system if the heat transfer problem between the plate and 
the flow is to be considered. In this work and consequently 0, is assumed to be 
held constant. 

If equations (78) and (79) are considered, the values of tin.' A x coefficients are 

(SI) 


(82) 
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Note that the relation A 2 = A\ y^/z/A holds if the shear 
glected, as in the case of linear piston theory. Therefore, 


stress effects are ne- 


(no shear) f 

•A^lpist 


2(l + o m ) + — (1 — a m ) y/0 


In practice 0 > 1, with the result that the previous ratio becomes always larger 
than approximately 4, irrespective of the value of a m . This suggests that piston 
theory estimates of panel stability are unconservative with respect to free molecule 
results when shear stress contributions are not considered. Actually, the same kind 
of relation between *4 2 and A\ is also true for the case of pure specular reflection 
(a m = 1). Again piston theory airloads are unconservative. 

Finally, note that the value of P x corresponding to the steady value of r is given 

by 

_ 1 ~ ot m 

^steady — Q / — 

yjns 


=» Pr 


jy r steady ^ 


1 


(85) 


3.2.2 Transition Regime: Particle Simulation Methods 

If intermolecular collisions are to be taken into account, there are basically two 
approaches: 1) try to solve the Boltzmann equation which is the fundamental 
equation of kinetic theory if only two-body collisions are important, or 2) model 
the gas flow at the molecular level^^’^^ - ^^. 

The Boltzmann equation is written in terms of the velocity distribution func- 
tion, with the independent variables being those of phase space (coordinates and 
velocities). One of the major difficulties in obtaining solutions for nontrivial gas 
flow problems in the transition regime in this case lies in the complicated structure 
of the so-called collision integral. One way to overcome this problem is to linearize 
the collision term for flows where the average speed and temperature exhibit little 
variations. Another approach consists of using alternative, simpler expressions for 
the collision integral, i.e., collision models. A somewhat different method consists 
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in expanding the velocity distribution function in terms of a small parameter. If 
this parameter is taken as the Knudsen number, such an expansion is known as 
the Chapman- Enkosg expansion. Both the Euler and Navier-Stokes equations of 
fluid mechanics can be derived in this way by retaining terms of up to zeroth and 
first order, respectively. The natural extension of the Chapman-Enskog theory is 
to carry the expansion to higher powers of Kn, which should provide formulations 
applicable further and further into the transition regime. The so-called Burnett 
equations correspond to carrying the Chapman-Enskog expansion to second-order 
accuracy, and there has been a major effort at Stanford University to calculate 
hypersonic Burnett shock solutions 

The alternative, as stated before, is to use some kind of particle simulation 
method, where the gas is modeled by some thousands or even millions of simu- 
lated molecules in a computer. The memory requirements are very large in this 
case since the velocity components, internal states, and position coordinates of the 
particles must be stored as the solution is advanced in time. On the other hand, 
such phenomena as dissociation and chemical reactions can be readily incorporated 
into the program as long as the physics of the processes can be described at the 
molecular level. The calculations are always unsteady, with the molecules being fol- 
lowed through representative collisions and boundary interactions. However, most 
applications up to the time of this writing have been restricted to obtaining steady 
state solutions as the large-time limit of the unsteady flow. One of the main reasons 
for this is the necessity to use some kind of averaging to obtain accurate macro- 
scopic properties such as pressure and density. If the solution is assumed to lie 
approaching a steady state, time averaging can be used. A cell network is necessary 
only in physical space, and then mainly to facilitate the sampling of properties and 
the choice of collision pairs. Usually the algorithms are described in such a way 
that the computation time is directly proportional to the number of molecules and 
cells in the simulated flow. This means that solutions are also very computationally 
intensive because of the large number of simulated particles. 

The most widely used particle simulation method is the direct simulation Monte 
Carlo (DSMC) method, as introduced by G. A. Bird^^l i n the 1960s. This method 
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has been constinuously advanced to the state where it is reliable, quite accurate and 
can handle complex problems [11 'l In principle, three-dimensional problems are 
not conceptually different from two- or one-dimensional cases since the collisions 
are always treated as a three-dimensional phenomenon. However, the CPU time 
and memory requirements for a typical three-dimensional application may be up to 
two orders of magnitude larger than those for a two-dimensional problem. One of 
the major drawbacks of the DSMC method is that the collision algorithm employed 
does not allow one to make effective use of supercomputers having a vector and/or 
parallel architecture. This can be better realized if it is known that an efficiently 
vectorized code is over an order of magnitude faster on a vector machine than 
the same code when limited to scalar execution. Therefore, the development of a 
selection rule governing collisions which is compatible with parallel decomposition 
would greatly improve the performance of such a method. This was accomplished 
by Baganoff and McDonald [112] , leading to the theses of McDonald ( 113 J, applicable 
to the Cray vector computer architecture (initially the Cray-2 arid afterwards the 
Cray-YMP), and Dagum [114] , who used a 32k processor Connection Machine. 



Chapter 4 
RESULTS 


The general solution for the system of equations (43) in free molecule flow is ob- 
tained following the procedure described by Dowell ^9). Thus, all the parameters 
(A, Moo, a m) 7, th, P z , R x , o/k t , the kind of curvature and H/h ) are speci- 
fied, with the q n being determined as functions of time T by numerical integration. 
The numerical integrator used is an adaptive-stepsize, fifth-order Runge-Kutta algo- 
rithm (-©^1, monitoring of local truncation error to ensure accuracy. Actually, 
the code is implemented in such a way to allow the user to specify a maximum value 
for the time step, with a default value of 9.95 x 10 -3 in terms of the nondimensional 
time T. No results are going to be shown, but several studies were made regarding 
the maximum size of the time step. If the default value is used, the actual time 
step usually remains close to 8 X 10 -4 . For a maximum value of the order of 10 -5 
or smaller, there is basically no change in the time step, while the final solutions 
remain the same. Therefore, it may be assumed that there are no mathematical 
instabilities associated with the numerical integration. Usually the interesting kinds 
of solutions which are sought are limit cycles, with the value of A being varied for 
a given set of parameters until such solutions are found. 

Some aspects of the problem can be discussed before any actual numerical re- 
sults are shown. First, it is useful to take a closer look at the integro-diflerential 
equation (31), especially the terms involving p*. Now recall that this distributed 
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load has a contribution proportional to W^. Then the term 

(j£ + p.) w, t 

leads to a contribution of the form (W,£ ) 2 , which is similar to the one present in the 
second-order piston theory formulation, equation (50). As explained before, this 
produces an overpressure, tending to push the panel into the cavity, i.e., the limit- 
cycle motion may be expected to be asymmetric with respect to its maximum and 
minimum amplitudes. This having been said, one may conclude that the minimum 
amplitude (in absolute value) is going to be larger than the maximum amplitude in 
the case of a flat panel. As will be seen, this is not the case, most probably because 
of the integral terms involving p £ which multiply As a matter of fact, these 

terms can be written as 

f W(£,T)d£ — W(£ y T) + terms involving W& * 

Jo 

meaning that the total membrane stress for a given deflection does not equal the 
one for the same deflection with a sign change. This further indicates an asymmetry 
in the limit-cycle motion. 

Second of all, an analysis of the contribution of the airloads to the coefficient a n 
of the generic equation (43) can give some hints about the variation of the frequency 
of flutter (or the main harmonic of a limit-cycle solution). Looking at the definition 
of a n in Appendix A, one can see that the presence of the term A$/2 > 0 indicates 
that the frequency of flutter is going to be larger than the one corresponding to 
the same problem when shear stress effects are neglected. If the influence of initial 
curvature is considered, an increase in the initial amplitude of the panel (H/h >0) 
may be expected to lead to an increase in the main harmonic in the case of sinusoidal 
curvature. The behavior should be the opposite for a constant curvature (parabolic) 
panel. 

Thirdly, suppose that the panel is fixed horizontally at both ends. It can be 
shown that, in the case of a uniform distributed longitudinal load like the one 
corresponding to P r , half of the panel is under tension while the other half is under 
compression. Then a large enough load may cause the panel to buckle. When 
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thi s happens, the buckling mode shape has a considerable contribution from the 
second mode, in contrast to classical buckling. In the latter case the buckling 
mode shape corresponds to the first mode solely. It should be recalled that only 
a simply-supported panel is being considered here. If a six-mode representation is 
used P x buck = 83.16, and the ratio between the contributions of the second and first 
modes to the buckling mode shape equals 0.2615. 

Finally, regarding a general solution q n , one can note from equation (43) that 
— q n is not a solution. This should be borne in mind when considering stable static 
deformation solutions, such as a buckled solution. Mathematically speaking this 
effect is associated with “quadratic” terms like a” q r q n or b* n * m qf n . By looking at 
the definition of the coefficients of equation (43) in Appendix A, it is clear that 
all such quadratic terms appear because of the unsteady shear stress and/or initial 
curvature. 

4.1 Nominal Configuration 

Regarding the numerical results to be presented in this chapter, the nominal con- 
figuration is specified by the following parameter values: 

M 0 0 = 25 , fi = 2.37 x 10- 9 , 

0 = 3.5 , th = 0.005 , 

a m = o , 7 = 1.4 , 

a h = 1- , «Jfc 2 = 1 , 

= 0, P x = 0, 

H/h = 0 , v = 0.3 . 

The environmental parameters correspond to a standard altitude of 110 km [110j , 
while the material is taken to be Rene 41. The value of 0 = 3.5 is associated with 
a panel temperature of approximately 840°K ~ 1050 °F. It is worth pointing out 
that these are realistic values corresponding to a point in the National Aerospace 
Plane (NASP) flight envelope [120,121] w hich approaches orbiting motion. 
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The first interesting observation when solving the problem is associated with 
the extremely small value for /i due to the low density. In this case, as mentioned 
in Reference [49], it is not possible to achieve a limit-cycle solution through time 
integration. Whatever motion is excited by the specified initial conditions, it con- 
tinues “forever” in the computation [Figure 4a)]. In order to achieve a limit-cycle 
solution, and thus define the flutter boundary, it was decided to introduce struc- 
tural damping of the viscous type into the formulation. This was accomplished by 
following the rule discussed in Section 2.3, in which As may be noted 

by comparing Figures 4a) and 4b), this indeed introduces enough damping into the 
system to allow the limit cycle to be achieved. The final limit cycle amplitude has 
clearly not been reached at T ~ 8, but the solution is changing toward it. The 
influence of the magnitude of this kind of damping on the flutter boundary is going 
to be discussed later. 

Next, it is also interesting to analyze the shear stress effects in the panel flutter 
phenomenon. Figure 5 shows the bifurcation diagram (limit-cycle amplitude versus 
dynamic pressure) for three different free-molecule models: FM p indicates that 

only the pressure is considered, FM US that unsteady shear stress effects are also 
included, while FM corresponds to the complete free-molecule model. One point 
to be mentioned here is that all values of the panel deflection to be shown in this 
chapter are associated with the coordinate £ = 0-75. This is appioximatelj the 
location where the panel assumes its maximum deflection amplitude in a given 
cycle. It should be said that a six mode representation is being used here; further 
comments on the number of modes and convergence of the solution are given later 
in this section. It can be seen from Figure 5a) that the introduction of shear stress 
effects stabilizes the panel. There is not a large difference if the steady part of 
the shear stress is considered or not, although its inclusion slightly destabilizes the 
system for small values of the limit-cycle amplitude. It is not shown here, but this 
is not true for an incompletely-converged two-mode solution, since the cuive for the 
FM (complete) model is significantly shifted to the right with respect to that for the 
FM US model. Another point to emphasize is that the limit-cycle peak amplitudes 
are slightly asymmetric due to the unsteady longitudinal load p^, as expected. This 
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a) No structural damping. 



b) Structural damping included: £1 = 0.01. 


Figure 4: Time history of the first generalized coordinate q x : nominal configuration, 
six-mode solution, A = 312.3 and ?i(0) = 0.01 [six mode representation, 
with all zero initial conditions except for ?i(0)]. 
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effect becomes more noticeable as the amplitude increases. In general the negative 
peak (in absolute value) is smaller than the positive peak. 

The variation of the limit-cycle main harmonic is shown in Figure 5b) for the 
different free-molecule models. Actually, all q n are essentially simple harmonic in 
the present case, although this may not be true in general. It should be noted that 
the frequencies used are made dimensionless in the same way as time. The usual 
natural frequencies for the linear panel correspond to the relation f n = n 2 7r/2. As 
expected, the inclusion of shear stress effects causes the frequencies to increase. 
A curious point is that the curves corresponding to the models FM P and FM US 
seem to be parallel, with the steady longitudinal load causing a sharper increase in 
frequency as the limit-cycle amplitude increases. In dimensional terms, the values 
of the frequencies in Hz are given by 

f 

For Rene 41 E a 1.8 x 10 n N/m 2 and p p « 8.2 x 10 3 kg/m 3 , such that for the 
nominal configuration f 7 f j a. Assuming a panel 60 cm long, this means that a 
value of / = 6 corresponds to a frequency of 70 Hz. 

Once it is understood how the inclusion of the aerodynamic shear stress affects 
the phenomenon, all remaining calculations were obtained using the most realistic 
FM or complete free-molecule model. The natural next step consists of a con- 
vergence study, illustrated in Figure 6 through the use of a bifurcation diagram. 
The numbers close to the solid curves indicate the numbers of modes used in each 
solution. It is clearly seen that six modes will ensure converged solutions, which 
compare very well with solutions obtained using up to twelve modes. The linear 
flutter parameter A cr , corresponding to a zero limit-cycle amplitude, can be de- 
termined with a relative error of approximately 15% when using only two modes, 
but this error increases for a given limit-cycle amplitude. For the six-mode solu- 
tion A cr = 312.296 and f cr — 5.6140. Other examples have been carried out which 
further verify these conclusions. 
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a) Limit-cycle amplitude. 



\h>p~k 


b) Limit-cycle main harmonic. 

Figure 5: Variation of limit cycle characteristics for different free-molecule models: 
nominal configuration and six-mode solution (see text for label defini- 
tions). 
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Figure 6: Bifurcation diagram: convergence study for nominal configuration. 

Figure 7 shows phase portraits for the two- and six-mode solutions, respectively. 
The numbers close to each curve indicate the generalized coordinates. The positive 
peak amplitudes in both cases are similar, and one can see how the higher modes 
influence the results. The arrows in the figures indicate the direction of motion, as 
well as points corresponding to a given time instant. From this last information, 
it is clear that even modes are practically in phase with each other, as well as odd 
modes. On the other hand, even modes are essentially out of phase relative to odd 
ones. Finally, it is worth emphasizing the initially unexpected result that the second 
mode is the dominant one, instead of the first mode. This is believed to occur due 
to the buckling of the panel under the uniform distributed load P x . In fact, if the 
value of Pxbuck given in the introduction of this chapter is used, the buckling load 
in the present case corresponds to a value of A « 247. This means that the panel 
is inside the buckling region when it starts fluttering. Although it is not shown 
because it is physically unrealistic, a phase portrait corresponding to Figure 7 but 
including only unsteady shear stress effects leads to a dominant first mode. From 

now on, all results to be presented correspond to a six-mode solution, unless stated 
otherwise. 
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a) two-mode solution: A = 350 and ( — ) = 0.629. 

\h J peak 



a) six-mode solution: A = 450 and ( — ) = 0.677. 

V h / peak 


Figure 7: Phase portrait of generalized coordinates: convergence study for nominal 
configuration (numbers identify the trajectories of different degrees of free- 
dom). 
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Figure 8: Panel displacement at the points in time of maximum deflections: nominal 
configuration, six-mode solution and A = 450. 

The panel motion is illustrated in Figure 8 at two different times for the case 
of Figure 7b), that is, A = 450. The points in time correspond to the maximum 
deflections (positive and negative) of the limit cycle. The motion resembles clas- 
sical panel flutter, Reference [49], with the distinguishing difference that the point 
of zero deflection is closer to the center of the panel. This is obviously due to the 
contribution of the second mode, which also moves the point of maximum deflec- 
tion farther back (£ ~ 0.8). It should be emphasized that, due to the amplitude 
asymmetry, the points of maximum deflection along the panel for the positive and 
negative peaks are slightly different. 


4.2 Initial Conditions: Nonuniqueness 

Before any further analysis of the phenomenon, it is necessary to say something 
about the dependence of the solutions on the initial conditions. If <7i(0) ^ 0, with 
all the other initial conditions set to zero, this dependency can be partially seen by 
analyzing the time histories for different values of <7i(0), Figure 9. Note that the 
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a) 9l (0) = 0.0001. 



b) 9,(0) = 0.01. 


Figure 9: Time histories of the first generalized coordinate </,: nominal configura- 
tion, A = 312.3 and zero structural damping. 
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solutions were computed with zero structural damping in order to permit, a bettor 
visualization. Also, A = 312.3 corresponding to a value just above the linear flutter 
boundary. In spite of the quite different nature of the time histories, it can be seen 
that some kind of modulation exists, especially for the cases of ryi(O) = 0.0001 and 
0.01. If the frequency content of these curves is analyzed, one finds that the power 
spectra of the generalized coordinates have very distinct peaks for two different 
values of frequency, which get closer together as q } (0) decreases. The first peak is 
generally greater than the second one, with the difference between them increasing 
for larger values of q x { 0). Actually, for the cases of <y,(0) = 0.0001 and 1 then- is 
essentially just one peak in the power spectra of q x if they are not plotted using a 
logarithm scale. For </,(0) = 0.0001 this happens because the peaks are too close 
together, corresponding to a nondimensional frequency of 5.59G (A/ = 0.09S). For 
<7i(0) = 1 the second peak is much smaller (two orders of magnitude) than the 
first one at / = 4.958. The dependency on the initial conditions is revealed even 
better, however, if the corresponding phase portraits are examined (Figure 10). 
It is important to point out that these results get much more complex as A is 
further increased, with a FFT analysis of the motion showing an increasingly broad 
spectrum of frequencies. 

The next step consists of ascertaining the influence of structural damping on 
the flutter boundary, as well as on the limit-cycle solutions. Regarding the latter, 
once structural damping is introduced in the system solutions for the various initial 
conditions (with the dynamic pressure parameter being hold constant), all converge 
to the same limit cycle, except for the case </,( 0) = 1. Figure 11 shows the' phase 
portraits of the generalized coordinates for the two kinds of limit cycles obtained. 
The different values of A correspond to solutions that yield approximately the same 
value of the maximum positive deflection at £ = 0.75. Again the numbers close to 
the curves indicate the generalized coordinates. Solution type I corresponds to the 
basic solution discussed in the last section. Solution type II represents a much more 
complex limit cycle, with the panel attaining much larger deflections: a maximum 
positive deflection at £ = 0.75 of 0.3100 versus 0.0024 for the other solution when 
A = 312.3. Also, this second kind of solution has a frequency content richer than the 
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Figure 10 : Phase portraits of the first generalized coordinate <7,: nominal configu- 
ration, A = 312.3 and zero structural damping. 
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Figure 10: (cont.) Phase portraits of the first generalized coordinate q x : 
configuration, A = 312.3 and zero structural damping. 


nominal 




AFTER 4. RESULTS 



nominal configuration 
of different degrees of 



4.3. EFFECTS OF AERODYNAMIC PARAMETERS 


55 


first one. Interesting to note is that the main harmonics are at 5.614 and 5.003 for 
solution types I and II, respectively. They compare very nicely with the frequencies 
obtained for the cases gi( 0 ) = 0.0001 and 1 when the structural damping was taken 
as zero. For the rest of the analysis in this chapter it is going to be assumed that the 
type-I limit cycle is the significant solution. This decision was based on two facts: 
1 ) <?i( 0 ) = 1 is quite a large “disturbance”, and 2 ) the type-II limit-cycle solution is 
not “robust” in the sense that such a solution cannot always be obtained for a given 
A if the initial conditions are the ones corresponding to the same kind of solution 
for a different value of the dynamic pressure parameter in the neighborhood of A. 
The type-I limit cycle is robust in this sense. 

Going back to the question of how the flutter boundary is affected by the mag- 
nitude of the structural damping, Figure 12 shows that this influence is minimal for 
values of < 0.01. This is interesting because physically one would expect the ac- 
tual structural damping to be of order 0.001 or smaller. A value of G = 0.01 was 
used in all the results to be shown in this chapter in order to make the numerical 
calculations converge faster. It is worth pointing out that larger values of ^ lead to 
larger values of A cr , as discussed in Section 2.3. For example, A cr ^ = 0 .ooo 5 = 312.286 
as compared to the value of 312.296 given previously. 


4.3 Effects of Aerodynamic Parameters 

There are basically three aerodynamic parameters in the present problem: 1 ) the 
Mach number M^, 2) the temperature ratio 0, and 3) the momentum accommoda- 
tion coefficient a m . The mass ratio /i could also have been included, but in the case 
of free-molecule flow it is too small to be of significance. The Mach number and 
mass ratio only appear as the combination n/M when piston theory is used. This 
is no longer true for free molecule flows, with the coefficient M 3 being proportional 
to the Mach number, equation (83). 
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Figure 13: Influence of the temperature ratio 0 on the linear flutter parameter A cr , 
and comparison with piston theory results. 


Comparison with piston theory results are shown in Figure 13 for the linear 
flutter parameter. The figure also shows how the free-molecule A cr varies with 
the temperature ratio for the two limiting cases of a m = 0 and 1. Note that in 
the latter case, when all the molecules are reflected specularly from the panel, the 
flutter condition is independent of 0. It is clear that the comparison between piston 
theory results and those obtained using the present rarefied formulation is highly 
dependent on the values of 0 and a m specified. The fact that A cr decreases with 0 
for a m < 1 should have been expected, since an increase in the temperature ratio 
implies an increase in the pressure load, as indicated by the expressions for A\ and 
A 2 , equations (81) and (82). 

The effect on the linear flutter parameter of varying a m is shown in Figure 14. 
The first thing that strikes one is that the curves are nearly straight lines for high 
enough values of the temperature ratio. This means that it may sufficiently accurate 
during a design process to only compute two points in the curve, and use them to 
estimate A cr for different values of a m . It is suggested in the literature that values 
in the vicinity of 0.2 lead to better agreement with experimental results. Another 
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Figure 14: Influence of the momentum accomodation coefficient a m on the linear 
flutter parameter X CT . 


point is that the slopes of the curves in Figure 14 depend on the value of 0 which 
is assumed. In fact, this could have been realized from Figure 13 since the free- 
molecule curves shown there represent limiting cases; that is, given a value of 0, A cr 
for an arbitrary value of a m must lie between the curves for a m = 0 and a m = 1. 
This gives rise to a curious observation: there is a value of the temperature ratio 
(0 « 5.6) for which A cr is independent of a m . 

Finally, Figure 15 gives the variation of the linear flutter parameter with Mach 
number. It is seen that A cr monotonically decreases with decreasing M 0 0 . If the 
effect of the terms associated with *4 2 and A 4 is neglected, which is a reasonable 
assumption due to the small value of the mass ratio, this trend makes sense in view 
of the fact that A 3 decreases with a decrease of the Mach number, and that the 
shear stress stabilizes the system. In the case of linear piston theory, A cr increases 
when Moo decreases for a fixed value of the mass ratio, that is, the trend is the 
opposite one. 
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Figure 15: Influence of the Mach number M <*> on the linear flutter parameter A cr . 

4.4 Effects of Structural Parameters 

Three set of parameters are included in this category: 1) the “spring support param- 
eters” a*,, 2) the applied compression parameter R x , and 3) the pressure differential 
parameter P z . 

In the first case, it is assumed that the panel is supported at leading and trailing 
edges by springs oriented longitudinally, such as to idealize flexibility in adjacent 
panels or the underlying structure. In classical panel flutter the presence of such 
springs alters the limit-cycle amplitude but not the linear flutter parameter. That is, 
only the supercritical characteristics of the panel are modified. In the present case 
the behavior is quite different because of the distributed tangential loads, especially 
P x . Under such a uniform distributed load, and for both ends fixed in the x direction, 
half of the panel is under tension (front), while the other half is under compression 
(back). When the trailing edge is supposed to have a spring attached to it (a kl = 1, 
n k 2 < 1), the smaller the spring constant, the larger the tension region along the 
panel. This appreciably stabilizes the system, as can be seen from Figure 1G. On the 
other hand, when only the leading edge spring exists («*., < 1, a k2 = 1) the behavior 
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Figure 16: Influence of the spring support parameters a k , on the linear flutter 
parameter A cr . 

changes. Now the smaller the spring constant, the larger the compression legion, 
and the system is destabilized. However, the rate of decrease of A cr is much smaller 
than the rate of increase of A cr in the first case. Finally, it is shown in Appendix 
B that the value of the linear flutter parameter for an arbitrary combination a kl , 
a k2 can be found by considering the appropriate case where only the front spring 
or the back spring exists, depending whether a ki > a k2 or not. 

Before commenting on the influences of other parameters, it is interesting to 
analyze the behaviors of the system when the a k , are varied for different values of 
the temperature ratio 0. In order to do this it is useful to define two new variables. 

A e = t ~ An ° — + 1 and A e = T“~ - ( S7 ) 

A n0 m ,0=3.5 ^nom.0 

where A nom 0 corresponds to the value of A cr for the nominal configuration but a 
given value of 0. Note that the second variable corresponds to a normalization of 
the dynamic pressure parameter for each different temperature ratio. On the other 
hand, the first variable is nothing more than a translation in the origin of the A 
scale, plus a normalization with respect to the nominal configuration. This last 
normalization is done only to make it possible to plot the flutter boundaries using 
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the same scale for both A© and A@. One should realize that A^ =3 5 and A© =3 .s are 
the same variable. Figure 17 shows how the linear flutter parameter varies with 
the spring support parameters using these variables as the ordinate. Three values 
of the temperature ratio were investigated, 0 = 1.5, 3.5, and 6.0. It is very clear 
that the use of \‘ e brings the flutter boundaries closer together, especially the lower 
curves corresponding to a kl < 1, <*k 2 = 1. Also, note that the flutter boundaries are 
steeper for smaller values of the temperature ratio. 

The behavior of the system in the presence of applied compressive loads, R x < 0, 
is illustrated in Figure 18. It resembles very much the kind of results obtained in 
classical panel flutter, with a distinguishing difference associated with the buckling 
boundary. As mentioned before, for a given set of parameters in the case of a static 
deflection, if q n is a solution, - q n is not. Computations show that the positive-static 
solution corresponding to a given A is smaller than the negative one (in absolute 
value). Also, as - R x increases beyond n 2 this difference increases, with the static- 
positive solution approaching zero earlier than the corresponding negative solution. 
Both these effects are illustrated in Figure 19. Then for large enough values of 
- R x the solution becomes noticeably “unstable” for values of the negative static 
deflection close to zero (Figure 20). The meaning for “unstable” here is that the 
solution seems to be converging with time to a finite value, but it then goes to zero 
in a relatively short interval. The interpretation of this result is that the system 
tries to approach a negative deflection when the corresponding positive solution 
has already become zero. Some disturbance then makes it “jump” from the former 
to the latter. It should be mentioned that the behavior of the first generalized 
coordinate (Figure 20) is typical of the entire motion. 

Another interesting feature related to the action of compressive loads is illus- 
trated in Figure 21, which is representative of points close to the intersection be- 
tween the buckling and flutter boundaries, R x « - 33. The case shown in the figure 
corresponds to a value of R x slightly smaller than the one associated with the inter- 
section between boundaries. Again the solution seems to be converging to a finite 
value, but it then begins to vary rapidly. In this case it eventually would converge 
to zero if the calculation were allowed to run for a longer time. However, if values 
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a) Shifted variable. 



b) Normalized variable. 


Figure 17. Influence of the temperature ratio 0 in the flutter boundaries when the 
spring support parameters a*, are varied. 
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Figure 18: Influence of the applied compression parameter R x on the linear flutter 
parameter A cr . 


of - R x larger than the one corresponding to the intersection between boundaries 
are used, a similar behavior can be noticed with the further complication that the 
motion does not seem to subside. This failure to converge to rest can also happen 
in the limit-cycle solution region, which was not the case for smaller values of - R x - 
As discussed by Dowell and Ilgamovl 59 ’ there is the possibility that this motion is 

an example of chaos. 

The influence of the temperature ratio on the flutter boundaries is again inves- 
tigated in Figure 22. As in the previous case (Figure 17) the stability boundaries 
come close together when the variable \* Q is used to plot the curves. In fact, the sta- 
bility boundaries for 0 = 3.5 and 6.0 practically coincide. Another interesting point 
is that the slope of the limit-cycle boundary, which is essentially linear, increases 
as the temperature ratio increases in this case [Figure 22b)]. On the other hand, 
this same slope decreases as 0 increases when A is used [Figure 22a)]. Basically 
this means that one obtains conservative results if A cr for any specified temperature 
ratio is estimated using the value of A* r 0 for the largest temperature ratio being 
considered. In other words, the value of the linear flutter parameter for a given 0 



CHAPTER 4. RESULTS 



50 100 150 200 250 


a) R x = — 7n 2 /6. 



b) R x = — 87t 2 /3. 


Figure 19: Influence of the applied compression parameter R x in the amplitude of 
the solution. 
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Figure 20: Time history of the first generalized coordinate q\ for R x = — Sir 2 / 3 and 
A - 92.656. 
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Figure 21: Time history of the first generalized coordinate q x for R x — — 107r 2 /3 
and A = 93.704. 
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obtained using as an estimate A* e for a larger temperature ratio is smaller than 
the correct value. 


Pressure differentials are specified by the parameter P z . A positive value corre- 
sponds to a uniform static pressure difference between the panel sides in the positive 
z direction, i.e., excess pressure in the interior. It can be seen from Figure 23 that 
the flutter boundaries again approach each other if they are plotted using Aq, i.e., 
their curvature are more nearly the same. As in the case of applied compressive 
loads, a large temperature ratio yields smaller values of A* r 0 when compared with 
the same results for a smaller temperature ratio. A new feature of the free-molecule 
case compared to classical panel flutter is that the flutter boundary is no longer 
symmetrical about its location for P z = 0. An upward pressure ( P z = 0) causes 
a somewhat lower increase in the linear flutter parameter than a downward pres- 
sure. Finally, it should be recalled that in the presence of a pressure differential 
the limit-cycle solution is “superimposed” upon a nonzero static deflection. At the 
linear flutter condition Figure 24 shows tha.t the deflection increases in absolute 
value with increasing values of P z . It can also be seen that the static deflection is 
essentially antisymmetric with respect to the line P z = 0. 


Finally, Figure 25 shows how the amplitude of the limit-cycle solution varies 
for different values of the temperature ratio. In order to understand the plots in 
the figure, it should be emphasized that none of the curves have the same value of 
A C r,e or A* r @, as seems to be implied by the plots. This feature was arrived at by 
appropriately translating the origins of the horizontal axes of the diagrams for the 
different values of 0 in order to obtain a better visualization. This accounts for the 
fact that there is no numbering in those axes. It can be seen from Figure 25 that 
the limit-cycle amplitude increases faster with Ae the higher the temperature ratio 
is, the difference in curvature being easily verified. However, if A @ is used instead 
the curves are much closer together, especially for larger values of 0. The same 
kinds of results are obtained if values of P z / 0 and ag ^ 1 are used. 
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-R z 

b) Normalized variable. 

Figure 22: Influence of the temperature ratio 0 in the stability boundaries when 
the applied compression parameter R x is varied. 




^r,e 0.8 



Figure 23: Influence of the temperature ratio 0 in the flutter boundaries when the 
pressure differential parameter P z is varied. 
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Figure 24: Influence of the pressure differential parameter P z on the static deflection 
at the linear flutter condition. 

4.5 Initial Curvature Effects 

When incorporating the effects of initial curvature into the code, a comparison with 
the results by Dowell ^2) was performed to ensure the correct implementation of the 
coefficients in the general equation of motion (43). As a direct result of this effort, 
it was established that the problem becomes very sensitive to the initial conditions 
as the maximum height of the panel, //, increases. Furthermore, the motion also 
tends to be what Dowell called a “heavily modulated oscillating function” for large 
values of the panel height. Most probably these motions constitute yet another 
example of chaos. Because of this fact, the present study was limited to relatively 
small values of the parameter H/h. 

Only results for two kinds of initial curvature are analyzed here: 1) constant 
curvature, and 2) half-sinusoidal curvature (u> 0 = H sin 7 rf). In these cases the 
undeformed shapes of the panel are not markedly different from each othei. Ac- 
cordingly, the bifurcation diagrams for both cases when Hjh — 0.5 are similar, with 
values of A cr close to each other (Figure 26). It should be noted that in both cases 
the system is destabilized with respect to the flat panel problem. This is the same 
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Figure 26: Bifurcation diagrams for H/h — 0.5. 

trend obtained when linear piston theory is used. For small values of A the solution 
corresponds to a static deflection due to the steady airloads caused by the panel 
curvature. The deflection (from the initial curved shape) at the point f = 0.75 is 
initially positive but then becomes negative. The panel deformation is illustrated in 
Figures 27 and 28 for two different values of A in each case. It can be seen from part 
a) of the figures that a large portion of the panel is deflected downward (relative to 
the initial curved shape) even for small values of the dynamic pressure parameter, 
and that this deformation approaches a half-sine shape as A is increased. On the 
other hand, if the deformation is superimposed to the initial shape [part b) of the 
figures], it becomes clear that the panel starts fluttering when it is deformed to a 
relatively flat geometry. In these conditions the compressive loads acting on the 
panel are, most probably, considerable. 

When the maximum height is further increased to a value of H/h = 1.0, the 
behaviors of the system for the different cases of initial curvature are no longer very 
similar, with the constant curvature panel having a much smaller value for A cr (Fig- 
ure 29). The static-deformation behavior of the panel resembles the one discussed 
for H/h = 0.5. One complication now is that, once the panel starts fluttering, 





a) Perturbed deflection. 



£ 

b) Total deflection. 


Figure 27: Static deflection of parabolic panel for H/h = 0.5. 
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Figure 29: Bifurcation diagrams for H/h = 1.0. 


the solution becomes very sensitive to the initial conditions. Eventually there is a 
change in the kind of the limit cycle achieved, even if A is increased through very 
small increments. 

The behavior of the system for a value of H/h — 2.0 is similar to the previous 
case of H/h = 1.0, with a further separation between the values of A cr for the two 
initial curvature cases discussed. It is interesting to observe, however, that the linear 
flutter parameter has increased for both cases relative to the values corresponding 
to the problem when H/h = 1. 

In view of the foregoing results, any actual analysis in a design process should 
be done assuming a parabolic shape, since this would give smaller values for the 
linear flutter parameter, leading to conservative results. On the other hand, due to 
the sensitivity of the results to the initial curved shape, one may wonder whether 
there really is some meaning in doing this kind of analysis, especially because of 
manufacturing restrictions on the panel final shape. 
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Figure 30: Bifurcation diagrams for Hjh = 2.0. 


4.6 Particle Simulation Method Application 


In principle, the interaction between the structural code and the aerodynamic one 
can be implemented as in other recent aeroelastic applications which use Com- 
putational Fluid Dynamics (CFD). The only question that should deserve some 
investigation is the kind of numerical integrator to be used, since a Runge-Kutta- 
like algorithm can no longer be used. However, this is not expected to be a serious 
problem. Therefore, the major difficulty in applying a particle simulation method 
to the present system has to do with CPU time requirements, since these methods 
are very computationally intensive. 

As mentioned before in Section 3.2.2, the CPU time can be assumed to be 
directly proportional to the number of molecules and cells in the simulated flow. 
Actually, the number of cells can be directly written in terms of number of time 
steps. In the case of the schemes implemented by McDonald and Dagum^^] 
(Stanford Particle Simulation Method) a time step in the code is typically of the 
order of 1/3 of the cell size. Then, in order to obtain an estimate of the CPU time 
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requirements related to the present application it is useful to define a new time 
variable: 



Note that t* is a nondimensional time which measures how many panel chords have 
been traveled. All symbols have been defined previously and can be identified from 
their definitions in the Nomenclature. However, it should be recalled that T is the 
nondimensional time used during the time integration. 

A reasonable value of T to obtain a limit cycle solution, 7^>i, cannot be expected 
to be much less than 10. On the other hand, from the results previously presented 
in this chapter it can be seen that the smallest possible value of A for a limit-cycle 
solution is of order 100. Because of restrictions on the applicability of particle 
simulation methods, the Mach number should not be taken much smaller than 10. 
Assuming a value of A M <*> = 10 3 , one can obtain the total number of chords traveled 
in order to ensure a converged flutter solution for different values of as shown in 
Table 1. 

Table 1: Variation with \i of the number of chords traveled to achieve a limit-cycle 
solution: A M <» = 10 3 and 7^ 0 i = 10. 



CO 

1 

O 

r— t 

lO " 4 

io - 5 


E 

O 

1 

chords traveled 

10 4 

3.16 x 10 4 

10 5 

3.16 x 10 5 

10 6 


The results in Table 1 can be put in terms of CPU time requirements. From 
Dagum’s work, CPU time requirements are of 2.0 /isec/particle/timestep when one 
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uses the 32k processor Connection Machine, and 1.7 /xsec/particle/timestep on 
the Cray 2. It is reasonable to assume that this time has decreased to at most 
1.0 //sec /particle / timestep on the Cray YMP; possibly a Connection Machine with 
a larger number of processors would require similar times. The largest possible 
number of particles should be used in order to assure time accurate results. On the 
other hand, the number of cells along the panel should be large enough to avoid 
excessive numerical errors when obtaining the aerodynamic generalized forces. A 
reasonable grid could be one with 250 cells in the x direction (about 200 cells along 
the panel), and 50 cells in the y direction. These values reflect the considerations 
arrived at by ^Voronowicz in his applications of the Stanford Method. The 
position of the wing leading edge, corresponding to the point where the boundary 
layer starts developing, should be at least a few cells ahead of the actual panel 
leading edge. On the other hand, the position of the wing trailing edge should be 
specified in such a way to allow for a sufficient number of cells after it, because of 
the subsonic region in the boundary layer. The extent of the grid in the y direction 
is necessary to avoid spurious reflections of shock waves at the upper boundary, for 
example. These values give a total number of 12,500 cells, and it follows that the 
smallest reasonable number of particles is approximately 4 x 10 6 , which allows for 
approximately 300 particles/cell. In the case of about 100 cells along the panel, 
which is a marginally small number, the number of particles might be decreased to 
about 10 6 if all relative proportions are kept approximately constant. As mentioned 
before, a typical time step is of the order of 1/3 of the cell size. Mdth these con- 
siderations, Tables 2 and 3 give Cray YMP CPU time requirements (in days) for 
the cases of 100 and 200 cells along the panel. It is clear that only the applications 
associated with the 100 cells case, and larger values of //, might become feasible if a 
decrease of two order of magnitudes could be arranged. The values in these tables 
were calculated by the following formula: 

10 -6 

CPU time = — — — (no. of particles) x 

24 x 3600 

x (value from Table 1) (3 x no. cells along the panel) . (89) 
Note that (value from Table 2) (3 x no. cells along the panel) gives the number 



78 


CHAPTER 4. RESULTS 


Table 2: Variation with n of CPU time requirements (days) to achieve a limit-cycle 
solution: A = 10 3 , T ao i = 10, and 100 cells. 



10" 3 

1 

O 

l—H 

io- 5 

10" 6 

io- 7 

days 

34.7 

109.8 

347.2 

1098 

3472 


Table 3: Variation with [i of CPU time requirements (days) to achieve a limit-cycle 
solution: A M '<» = 10 3 , T bo i = 10, and 200 cells. 



IO" 3 

IO" 4 

l—i 

o 

1 

Cn 

10" 6 

io- 7 

days 

277.8 

878.4 

2778 

8784 

27778 


of time steps for T so \ = 10. 

The first observation to be made if one is to investigate the possibility of smaller 
CPU times is that the values in Tables 2 and 3 are proportional to yj A / /i . This 

is the parameter in the estimates which may be directly affected by the panel’s 
characteristics (material properties). However, since it is obtained from nondimen- 
sional parameters of the panel flutter problem, a change in the stiffness of the panel 
for example is not reflected in \J A M m / ^/, and consequently in the CPU time re- 
quirements. Naturally, the physical time t is greatly affected. On the other hand, a 
lighter material makes larger values of p more reasonable. As a limiting case, con- 
sider corkwood, which has a density of approximately 0.21 g/cm 3 . This corresponds 
to almost a two-order-of-magnitude decrease when compared with the value of 8.20 
g/cm 3 for Rene 41. If an altitude of 95 km is assumed instead of 110 km, which 
was used in the calculations shown in this chapter, a further decrease of two orders 
of magnitude is obtained. This makes a value of fi = 10 -4 “reasonable”, but there 
is little hope of managing to decrease the CPU time with a further decrease in p 
and still have some physical material. Since the value of is also “fixed” close 
to 10, only A remains to be varied, that is, for smaller values of CPU time A should 
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be decreased. Of all the nondimensional parameters involved in the panel flutter 
calculations, only four of them lead to a decrease of the linear flutter parameter 
for the flat panel problem. They are the Mach number, the temperature ratio, the 
spring parameter on the leading edge, and the applied compression parameter. By 
far the most strong influence is related to increasing 0. A complete study was done 
in order to find conditions which would lead to smaller value of A, but the minimum 
value obtained for A cr was approximately 35. 1 This is almost a decrease of one order 
of magnitude in A, but it is definitely inadequate since CPU times are proportional 
to vX 

The final conclusion of this brief investigation is that, until supercomputer speeds 
are increased by a factor of at least 100 it will be infeasible to carry out flutter 
calculations with particle-method aerodynamics. 


1 A value of the temperature ratio of 36 was used, which at the assumed altitude of 95 km means 
that the panel temperature is of the order of 6, 780 0 K « 6, 500 0 C « 11,700° F. This is quite high 
to say the least. 



Chapter 5 


FINAL REMARKS 


The main accomplishment of this work corresponds to the finding that aerodynamic 
shear effects are considerable in the rarefied (free-molecule) flow analyzed. Here the 
ratio between shear stress and pressure, for both steady and unsteady parts of 
the loading (up to first order), is of order of Mach number (see page 34). One 
consequence of including the aerodynamic shear loading is to stabilize the panel 
response, with an increase of the linear flutter parameter* of 53% relative to the 
pressure only value in the case of the nominal configuration specified in Section 
4.1. Using results given in Section 6.5 of the book by Anderson [122], j s p OSS ible 
to show that in the case of a flat plate r/p^ ~ M^/yfRe^. Here r is the value 
of the steady shear stress at a given point of the plate, p^ is the static pressure 
corresponding to the freestream, and Re x is the Reynolds number based on the 
coordinate x measured from the leading edge of the plate. Then, depending on 
the position of the panel relative to the leading edge of the wing and the value 
of the Mach number, the ratio r/p may be large enough to change the flutter 
characteristics of the panel significantly. 

A convergence study regarding the number of modes to be used in the Galerkin 
method showed that a six-mode solution leads to converged results. Up to 12 
modes were used in this research. The motion resembles continuum panel flutter, 

tThis parameter corresponds to the dynamic pressure parameter at the linear flutter boundary, 
that is, a limit cycle of zero amplitude. 
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with the distinguishing difference that the point of zero deflection (the node of 
the flutter mode) is closer to the center of the panel. The point of max, mum 
deflection is also slightly farther back (4 « 0.8). Both these changes are a d.rect 
consequence of the fact that the second mode is the dominant one, as opposed to 
the first mode in the case of classical panel flutter. It is believed that this happens 
because of the interaction between flutter and buckling in the presence of a uniformly 
distributed longitudinal load. Finally, it must be mentioned that the limit-cycle 
peak amplitudes are slightly asymmetric due to the unsteady longitudinal load. 

Some key findings about the effects of the several parameters considered are as 

follows: 


Temperature ratio 0 = Tp/T*,: this is the ratio of panel temperature to that 
of the freestream flow. The larger 0, the larger are the aerodynamic pressures. 
Accordingly, the linear flutter parameter decreases with increasing 0, all ot lici 
parameters being fixed. Their relationship may be quite nonlinear. 

“Momentum accommodation coefficient” o m : this number gives the fraction of 
molecules which reflects specularly upon collision with the panel. A value of 0 
implies that all collisions contribute to aerodynamic shear, with the molecules 
being reemitted in random or “diffuse” fashion according to the temperature of 
the panel. Calculations show that the variation of the linear flutter parameter 
with o m is close to linear if the temperature ratio is sufficiently large. In the 
case of the nominal configuration specified in Section 4.1, this comes about 
vpIiips larger than 0 — 3.5. 


. Mach number A#*,: the linear flutter parameter decreases monotonically with 
decreasing Mach number. This is contrary to the case of continuum flow 
if linear piston theory is used, where increasing the Mach number leads to 
smaller values of the linear flutter parameter. 


• Spring support parameters a ki : in this analysis it is assumed that the panel 
is supported at leading and trailing edges by springs oriented longitudinally, 
with spring constants k, and k 2 respectively. In classical panel flutter, such 
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springs would affect only the supercritical limit-cycle amplitude. However, 
their presence also changes the linear flutter parameter in free- molecule flow 
because of the influence of distributed longitudinal loads representing the 
aerodynamic shear. The panel is strongly stabilized if only the spring at the 
trailing edge is considered, the effect being the opposite if only the leading 
edge spring constant is varied. However, the rate of decrease of the linear 
flutter parameter is not as large as the rate of increase in the former case. 

• Applied compression parameter R x : this quantity relates the compressive end- 
load to its value that gives rise to classical column buckling. It has been found 
that the relationship between the flutter boundary and the buckling boundary 
does not change much from that which is typical of continuum panel flutter. 
Naturally, numerical values are different from case to case. 

• Pressure differential parameter P z : this corresponds to a uniform static pres- 
sure difference between the panel’s sides. An upward pressure (i.e., excess 
pressure in the interior) is associated with P z > 0. Regardless of the direction 
of the force, there is an increase in the linear flutter parameter with P z . The 
new feature of the free-molecule case is that the stability-boundary location is 
no longer symmetrical about the line P z = 0, with an upward pressure causing 
3 somewhat lower change than a downward pressure. 

• Initial curvature: the amount of curvature is specified by the ratio between 
the panel s maximum initial height and its thickness, H/h. Both a sinusoidal 
and a constant (parabolic) curvature panel were considered. For small values 
of H/h the flutter solutions for the two different cases are very similar. How- 
ever, for larger values of H/h the linear flutter parameter for the two kinds of 
curvature can be significantly different. This effect is already noticeable for a 
value of H/h = 1 in the cases that were studied in Chapter 4. The parabolic 
panel leads to smaller values of the linear flutter parameter, and so this kind 
of initial shape should be used in a design process if one would prefer to be 
on the conservative side. This last remark is made in view of the fact that 
the two shapes are not markedly different, with a maximum difference in their 
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ordinates of about 10% relative to H. Accordingly, for a very shallow curved 
panel (small ratio between the panel’s thickness and length) it may be unre- 
alistic to try to distinguish between one and the other due to manufacturing 
restrictions. 

There is a substantial effect of the temperature ratio on the location of the sta- 
bility boundaries for parametric studies related to the parameters Ofc, , R x , and P z . 
Different locations can, however, be “collapsed” onto one another by using as ordi- 
nate an appropriately normalized dynamic pressure parameter, given by equation 
(87). This works better for higher values of the temperature ratio. In most cases the 
linear flutter parameter obtained using this normalized definition decreases with in- 
creasing 0. This implies that conservative results (smaller values) can be obtained 
using as an estimate the linear flutter (normalized) parameter corresponding to 
the highest temperature ratio to be considered, and then transfoiming it foi each 
specific case of 0. 

Finally, a brief investigation was conducted on the feasibility of using a particle 
simulation method to obtain the aerodynamic loads. This resulted in the conclu- 
sion that such calculations can only become possible with a two orders of magnitude 
increase above the present speeds of supercomputers. Such a study would consti- 
tute a natural extension of the foregoing work, filling the gap between continuum 
and free- molecule flow. It is, therefore, recommended to a future generation of 
aeroelasticians. 


Appendix A 

GALERKIN INTEGRALS AND 
COEFFICIENTS 


A.l Sine and Cosine Integrals 

Some of the necessary integrals in the evaluation of the Galerkin coefficients arc 
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A. 2 Galerkin Coefficients 


Assuming that the parameter p specifies the kind of initial deformation, with p = 0 
corresponding to the constant curvature case, and p > 0 to the pth sine curve, the 
coefficients specified in equation (43) are given by 
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Appendix B 

LONGITUDINAL SPRINGS 
EFFECTS 


The influence of the longitudinal springs on the linear flutter parameter is given by 
the term multiplying P x on the expression for a n in Appendix A, which is 


P = (l-2— + a k 
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Assume ot kx /oi k2 = Then ot k can be rewritten as 
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where C * = 1/C. T can then be expressed in terms of a kl or a k2 alone: 
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First thing to note is how T and a k vary when the value of C changes. Then 
1. C > 1 a kl > a k2 . Then: 0 < T < 1 and a k < a kl . 

f90l 
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2. C < 1 => Q kj > «*,• Then: -1 < T < 0 and at k < a k2 . 

Secondly, the value of A cr for an arbitrary combination a kl , <Xk 2 can be found by 
considering one of the limiting cases a ki = 1, at k2 < 1 or ct k2 = 1, o kl < 1, depending 
whether d^, ^ Qt k 2 or not. Here ( • ) indicates an equivalent problem with different 
values for a kt . 


1. a kl > a k2 => d/t, = 1 , = 1/C 

The value of C can be found by equating J- and J 7 , that is, 


Then 
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(94) 


It can also be shown that a k < djt, which implies that the limit-cycle ampli- 
tude for the pair a ki , a k2 is larger than the one corresponding to a ki , a k2 if A 
is held constant. 


2. d fcl < Ofc 2 =» dfc 2 = 1 , a kl — 1/C 

The value of d^ can be found as in the previous case, except that now T and 
T should be written in terms of a k2 . Then 


<*ki 


(*k 2 ~ 2 
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(95) 


Again it can be shown that the limit cycle amplitude for the pair o^ , Ojt 2 is 
larger than the one corresponding to d*.,, dfc 2 if A is held constant. 
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